[GRASS-user] v.rast.stats
Sandile Gumede
akasandile at gmail.com
Wed Jul 7 09:15:43 EDT 2010
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20100707/4d75b39b/attachment-0001.html
More information about the grass-user
mailing list