[GRASS-user] Re: v.rast.stats
Sandile Gumede
akasandile at gmail.com
Tue Jul 13 07:19:58 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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/grass-user/attachments/20100713/ad39b1c2/attachment-0001.html
More information about the grass-user
mailing list