[GRASS-user] v.rast.stats
Sandile Gumede
akasandile at gmail.com
Fri Jul 9 08:53:02 EDT 2010
Is there any one who can help me?
On Wed, Jul 7, 2010 at 3:15 PM, Sandile Gumede <akasandile at gmail.com> wrote:
> Hi,
>
> I've tried to sort out the issue of the region resolution. Now when I'm
> running v.rast.stats it says:
>
> ERROR: No categories found in raster map
>
> Here is the script I'm running thats giving me that error:
>
> #!/bin/sh
>
> #variable to customize:
> # path to GRASS software main directory
> GISBASE=/usr/lib/grass64
> # path to GRASS database
> GISDBASE=$HOME/grassdata/Cape_Town
>
> LOCATION_NAME=SRTMDEM
> MAPSET=PERMANENT
>
> # nothing to change below
> MAP=$1
> LOCATION=$2
>
>
> # generate temporal LOCATION:
> TEMPDIR=FLOODS
> mkdir -p $GISDBASE/$LOCATION_NAME/$MAPSET
>
> # save existing $HOME/.grassrc6
> if test -e $HOME/.grassrc6 ; then
> mv $HOME/.grassrc6 /tmp/$TEMPDIR.grassrc6
> fi
>
> echo "LOCATION_NAME: $LOCATION_NAME" > $HOME/.grassrc6
> echo "MAPSET:$MAPSET" >> $HOME/.grassrc6
> echo "DIGITIZER: none" >> $HOME/.grassrc6
> echo "GISDBASE: $GISDBASE" >> $HOME/.grassrc6
> export GISBASE=$GISBASE
>
> # Create a WIND file with minimal information and no projection:
> echo "proj: 3
> zone: 0
> north: 1
> south: 0
> east: 1
> west: 0
> cols: 1
> rows: 1
> e-w resol: 1
> n-s resol: 1
> top: 1
> bottom: 0
> cols3: 1
> rows3: 1
> depths: 1
> e-w resol3: 1
> n-s resol3: 1
> t-b resol: 1
> " > $GISDBASE/$LOCATION_NAME/$MAPSET/WIND
>
> # Copy WIND-file to DEFAULT_WIND:
> cp $GISDBASE/$LOCATION_NAME/$MAPSET/WIND \
> $GISDBASE/$LOCATION_NAME/$MAPSET/DEFAULT_WIND
>
>
> echo "name: Latitude-Longitude
> datum: wgs84
> towgs84: 0.000,0.000,0.000
> proj: ll
> ellps: wgs84
> "> $GISDBASE/$LOCATION_NAME/$MAPSET/PROJ_INFO
>
> echo "unit: degree
> ubits: degrees
> meters: 1.0
> "> $GISDBASE/$LOCATION_NAME/$MAPSET/PROJ_UNITS
>
>
>
> export PATH=$GISBASE/bin:$GISBASE/scripts:$PATH
> export LD_LIBRARY_PATH=$GISBASE/lib:$LD_LIBRARY_PATH
> export GIS_LOCK=$$
> export GISRC=$HOME/.grassrc6
>
>
> # this should print GRASS version used:
> g.version
> # other calculations go here....
>
> # import rainfall data set.
> cd /home/tgumede1/grassdata/Cape_Town
>
>
> # rainfall data set.
> r.in.gdal input=$HOME/grassdata/Cape_Town/TRMMLast1day.tif output=rainfall
>
>
> # DEM data set.
> r.in.gdal input=$HOME/grassdata/Cape_Town/Dem_CF.tif target=SRTMDEM
> output=dem
> g.region rast=dem -p
>
> # creating set of maps indicating flow acc, drainage dir, streams
> r.watershed --o elevation=dem at PERMANENT drainage=flow_direction
> basin=catch accumulation=acc threshold=1 memory=300 stream=str
>
> # convert catch raster to polygon vector
> r.to.vect in=catch at PERMANENT out=catchments feature=area
>
> g.region rast=rainfall -p
>
> # Calculate univariate statistics
> v.rast.stats -c vector=catchments at PERMANENT raster=rainfall at PERMANENTcolpre=precp
>
>
>
> On Wed, Jul 7, 2010 at 9:18 AM, Sandile Gumede <akasandile at gmail.com>wrote:
>
>> Hi
>>
>> Which module do I use to change the resolutions?
>>
>>
>> 2010/7/6 <micha at arava.co.il>
>>
>>> Hello Sandile:
>>> It seems you are importing two raster with vastly different resolutions.
>>> (I think we already came across this). See below...
>>>
>>> > Hi
>>> >
>>> > Below is a step-by-step of what I have done but I'm getting an error
>>> when
>>> > running v.rast.stats vector=catchments raster=rainfall layer=1
>>> > colprefix=area.
>>> >
>>> >
>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > r.in.gdal
>>> > in=/home/tgumede1/grassdata/Cape_Town/TRMMLast1day.tif out=rainfall
>>> >
>>> > Projection of input dataset and current location appear to match
>>> > 100%
>>> > r.in.gdal complete. Raster map <rainfall> created.
>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > g.region rast=rainfall -p
>>> > projection: 3 (Latitude-Longitude)
>>> > zone: 0
>>> > datum: wgs84
>>> > ellipsoid: wgs84
>>> > north: 33:30S
>>> > south: 33:45S
>>> > west: 18:15E
>>> > east: 19E
>>> > nsres: 0:15
>>> > ewres: 0:15
>>> > rows: 1
>>> > cols: 3
>>> > cells: 3
>>>
>>> Here, the rainfall data has a resolution of 0:15 = 15 minutes or 1/4
>>> degree. THat's approximately (at the equator) about 27 km. So *one*
>>> raster
>>> cell is 27 km X 27 km =~ 730 sq.km. Your region is covered by 3 cells, 1
>>> row by 3 columns. Not very helpful data!
>>>
>>> Next...
>>>
>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > r.in.gdal
>>> > in=/home/tgumede1/grassdata/Cape_Town/Dem_CF.tif out=dem target=SRTMDEM
>>> >
>>> >
>>> > Projection of input dataset and current location appear to match
>>> > 100%
>>> > r.in.gdal complete. Raster map <dem> created.
>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > g.region rast=dem -p
>>> > projection: 3 (Latitude-Longitude)
>>> > zone: 0
>>> > datum: wgs84
>>> > ellipsoid: wgs84
>>> > north: 33:40:46.499215S
>>> > south: 34:00:52.499215S
>>> > west: 18:17:55.500436E
>>> > east: 19:10:16.500436E
>>> > nsres: 0:00:03
>>> > ewres: 0:00:03
>>> > rows: 402
>>> > cols: 1047
>>> > cells: 420894
>>>
>>> Your DEM layer, on the other hand, is of resolution 3 arc seconds, or
>>> about 90 meters on a side. So each cell is 90 m. X 90 m = 8100 sq.m. =~
>>> 0.0081 sq.km.
>>>
>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > r.watershed elevation=dem
>>> > accumulation=acc drainage=direction basin=catch stream=str
>>> threshold=200
>>> >
>>>
>>> Here you chose a threshold of 200. That's 200 cells, so about 200 X
>>> 0.0081
>>> sq.km, or 1.6 sq km. As a result (see below, the output of r.to.vect)
>>> you
>>> are getting over 19,000 tiny little catchments. Are you sure that's what
>>> you want??
>>>
>>> Finally, you're trying to get raster values for 19,000 tiny vector areas
>>> where the raster (rainfall) is only 3 cells! You'll have 1000's of
>>> catchments with all the same values. And I guess that some of these
>>> catchments are extending outside of the three rainfall cells, and causing
>>> the NULL value error.
>>>
>>> In short: I think you'll need to match the resolution of the DEM to that
>>> of the rainfall data. If the rainfall is only as accurate as 1 data value
>>> per 730 sq.km. then you will be able to do vector-raster analyses only
>>> at
>>> that resolution = i.e. continent scale maps.
>>>
>>> HTH
>>> --
>>> Micha
>>>
>>> >
>>> > SECTION 1a (of 5): Initiating Memory.
>>> > SECTION 1b (of 5): Determining Offmap Flow.
>>> > 100%
>>> > SECTION 2: A * Search.
>>> > 100%
>>> > SECTION 3: Accumulating Surface Flow.
>>> > 100%
>>> > SECTION 4: Watershed determination.
>>> > 100%
>>> > SECTION 5: Closing Maps.
>>> > 100%
>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > r.to.vect -s in=catch out=catchments
>>> > feature=area
>>> > Extracting areas...
>>> > 100%
>>> > 100%
>>> > Building topology for vector map <catchments>...
>>> > Registering primitives...
>>> > 60653 primitives registered
>>> > 314051 vertices registered
>>> > Building areas...
>>> > 100%
>>> > 19885 areas built
>>> > 1 isles built
>>> > Attaching islands...
>>> > 100%
>>> > Attaching centroids...
>>> > 100%
>>> > Number of nodes: 40769
>>> > Number of primitives: 60653
>>> > Number of points: 0
>>> > Number of lines: 0
>>> > Number of boundaries: 40768
>>> > Number of centroids: 19885
>>> > Number of areas: 19885
>>> > Number of isles: 1
>>> > r.to.vect complete.
>>> >
>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > v.rast.stats vector=catchments
>>> > raster=rainfall layer=1 colprefix=precip
>>> >
>>> > DBMI-DBF driver error:
>>> > SQL parser error: syntax error, unexpected NULL_VALUE processing 'NULL'
>>> > in statement:
>>> > UPDATE catchments SET precip_min=-NULL WHERE cat=10163
>>> > Error in db_execute_immediate()
>>> >
>>> > ERROR: Error while executing: 'UPDATE catchments SET precip_min=-NULL
>>> > WHERE
>>> > cat=10163'
>>> >
>>> >
>>> >
>>> > Here is the output of gdalinfo TRMMLast1day.tif
>>> >
>>> > Origin = (18.250000000000000,-33.500000000000000)
>>> > Pixel Size = (0.250000000000000,-0.250000000000000)
>>> >
>>> > --------------------
>>> > coordinates--------------------------------------------
>>> > Corner Coordinates:
>>> > Upper Left ( 18.2500000, -33.5000000)
>>> > Lower Left ( 18.2500000, -33.7500000)
>>> > Upper Right ( 19.0000000, -33.5000000)
>>> > Lower Right ( 19.0000000, -33.7500000)
>>> > Center ( 18.6250000, -33.6250000)
>>> >
>>> >
>>> > Here is what I did to clip the region of interest.
>>> >
>>> > gdal_translate -a_srs EPSG:4326 -projwin 18.2987501 -33.6795831
>>> 19.1712501
>>> > -34.0141665 3B42RT.2010032900.1day.tif TRMMLast1day.tif
>>> >
>>> >
>>> > Is there something I have done wrong in these steps or there is
>>> something
>>> > wrong with my coordinates?
>>> >
>>> >
>>> >
>>> > 2010/7/6 <micha at arava.co.il>
>>> >
>>> >> > Hi
>>> >> >
>>> >> > Is it wrong to use -a_ullr option in gdal_translate to clip a small
>>> >> > portion
>>> >> > of the region from the big geotiff file region?
>>> >> >
>>> >>
>>> >> The option -a_ullr will change the georeference of the resulting file,
>>> >> so
>>> >> you could say it's "wrong" if you want to keep the original
>>> referencing.
>>> >> The way to clip a portion of the original and still maintain
>>> >> geo-referencing is with the -projwin option.
>>> >>
>>> >> --
>>> >> Micha
>>> >>
>>> >> > --
>>> >> > Kind Regards
>>> >> > TS Gumede
>>> >> > CSIR, Meraka Institute
>>> >> > 072 258 1650
>>> >> >
>>> >> > This mail was received via Mail-SeCure System.
>>> >> >
>>> >> >
>>> >> >
>>> >>
>>> >>
>>> >
>>> >
>>> > --
>>> > Kind Regards
>>> > TS Gumede
>>> > CSIR, Meraka Institute
>>> > 072 258 1650
>>> >
>>> > This mail was received via Mail-SeCure System.
>>> >
>>> >
>>> >
>>>
>>>
>>
>>
>> --
>> Kind Regards
>> TS Gumede
>> CSIR, Meraka Institute
>> 072 258 1650
>>
>>
>
>
> --
> Kind Regards
> TS Gumede
> CSIR, Meraka Institute
> 072 258 1650
>
>
--
Kind Regards
TS Gumede
CSIR, Meraka Institute
072 258 1650
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20100709/78502709/attachment.html
More information about the grass-user
mailing list