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