[GRASS-user] Re: v.rast.stats
micha at arava.co.il
micha at arava.co.il
Fri Jul 9 17:20:37 EDT 2010
> 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.
>
>
>
More information about the grass-user
mailing list