[GRASS-user] v.rast.stats
Sandile Gumede
akasandile at gmail.com
Tue Jul 13 07:20:43 EDT 2010
Hi,
Thanks for the suggestions. I have downloaded the GTOPO30 data from their
website, its a compressed file e020n40.tar.gz I guess these are the
coordinates for part of Africa. I have decompressed the file and inside
there are 8 files; .DEM, .DMW, .GIF, .HDR, .PRJ, .SCH, .STX, .SRC.
Is the .GIF file the one I'm supposed to project in order to use it in
GRASS?
> 2010/7/9 <micha at arava.co.il>
>
> > These are my data sets attached:
>> >
>> > 1. Dem_CF.tif is the DEM
>> > 2. TRMMLast1day.tif is the rainfall data
>> >
>> > On Fri, Jul 9, 2010 at 2:53 PM, Sandile Gumede <akasandile at gmail.com>
>> > wrote:
>> >
>> >> Is there any one who can help me?
>>
>> I'm away this week, so I can only add a quick response:
>>
>> I still see two "glaring" problems in your technique.
>> 1- You mentioned that you try to "sort out the resolution" between the two
>> data sources. One source, the rainfall data, is a raster of 27 km. pixel
>> size, while the other, the DEM is 90 m. pixel size. I don't believe that
>> you can resolve this gap. You can not compare data with such a gap in
>> resolution. The only thing possible with rainfall data at that resolution
>> is *continent scale* maps. You can get the gtopo30 DEM (1 km. resolution,
>> I think) for all of Africa for example, and compare that with rainfall
>> data for all of Africa.
>> 2- When you get to the step of watershed creation, you seem to be using
>> illogical threshold values. Below you set threshold=1. That means every
>> single cell in the dem will become a basin! That's certainly *not* what
>> you want. Typical threshold values will be in the 1000's, even 10's of
>> thousands.
>> --
>> Micha
>> >>
>> >>
>> >> 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
>> >>
>> >>
>> >
>> >
>> > --
>> > 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/20100713/b61ac8fa/attachment-0001.html
More information about the grass-user
mailing list