Hi,<br><br>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.<br>
<br>Is the .GIF file the one I'm supposed to project in order to use it in GRASS?<br>
<br><div class="gmail_quote">2010/7/9 <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;">
<div>> These are my data sets attached:<br>
><br>
> 1. Dem_CF.tif is the DEM<br>
> 2. TRMMLast1day.tif is the rainfall data<br>
><br>
> On Fri, Jul 9, 2010 at 2:53 PM, Sandile Gumede <<a href="mailto:akasandile@gmail.com" target="_blank">akasandile@gmail.com</a>><br>
> wrote:<br>
><br>
>> Is there any one who can help me?<br>
<br>
</div><div>I'm away this week, so I can only add a quick response:<br>
<br>
I still see two "glaring" problems in your technique.<br>
1- You mentioned that you try to "sort out the resolution" between the two<br>
data sources. One source, the rainfall data, is a raster of 27 km. pixel<br>
size, while the other, the DEM is 90 m. pixel size. I don't believe that<br>
you can resolve this gap. You can not compare data with such a gap in<br>
resolution. The only thing possible with rainfall data at that resolution<br>
is *continent scale* maps. You can get the gtopo30 DEM (1 km. resolution,<br>
I think) for all of Africa for example, and compare that with rainfall<br>
data for all of Africa.<br>
2- When you get to the step of watershed creation, you seem to be using<br>
illogical threshold values. Below you set threshold=1. That means every<br>
single cell in the dem will become a basin! That's certainly *not* what<br>
you want. Typical threshold values will be in the 1000's, even 10's of<br>
thousands.<br>
--<br>
Micha<br>
>><br>
>><br>
</div><div><div></div><div>>> On Wed, Jul 7, 2010 at 3:15 PM, Sandile Gumede<br>
>> <<a href="mailto:akasandile@gmail.com" target="_blank">akasandile@gmail.com</a>>wrote:<br>
>><br>
>>> Hi,<br>
>>><br>
>>> I've tried to sort out the issue of the region resolution. Now when I'm<br>
>>> 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<br>
>>> output=rainfall<br>
>>><br>
>>><br>
>>> # DEM data set.<br>
>>> r.in.gdal input=$HOME/grassdata/Cape_Town/Dem_CF.tif target=SRTMDEM<br>
>>> 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<br>
>>> 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<br>
</div></div>>>> raster=rainfall@PERMANENTcolpre=precp<br>
<div><div></div><div>>>><br>
>>><br>
>>><br>
>>> On Wed, Jul 7, 2010 at 9:18 AM, Sandile Gumede<br>
>>> <<a href="mailto:akasandile@gmail.com" target="_blank">akasandile@gmail.com</a>>wrote:<br>
>>><br>
>>>> Hi<br>
>>>><br>
>>>> Which module do I use to change the resolutions?<br>
>>>><br>
>>>><br>
>>>> 2010/7/6 <<a href="mailto:micha@arava.co.il" target="_blank">micha@arava.co.il</a>><br>
>>>><br>
>>>>> Hello Sandile:<br>
>>>>> It seems you are importing two raster with vastly different<br>
>>>>> resolutions.<br>
>>>>> (I think we already came across this). See below...<br>
>>>>><br>
>>>>> > Hi<br>
>>>>> ><br>
>>>>> > Below is a step-by-step of what I have done but I'm getting an<br>
>>>>> error<br>
>>>>> 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>
>>>>> 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*<br>
>>>>> 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<br>
>>>>> cells,<br>
>>>>> 1<br>
>>>>> row by 3 columns. Not very helpful data!<br>
>>>>><br>
>>>>> Next...<br>
>>>>><br>
>>>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > r.in.gdal<br>
>>>>> > in=/home/tgumede1/grassdata/Cape_Town/Dem_CF.tif out=dem<br>
>>>>> 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>
>>>>> 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>
>>>>> =~<br>
>>>>> 0.0081 <a href="http://sq.km" target="_blank">sq.km</a>.<br>
>>>>><br>
>>>>> > GRASS 6.4.0RC5+39438 (SRTMDEM):~ > r.watershed elevation=dem<br>
>>>>> > accumulation=acc drainage=direction basin=catch stream=str<br>
>>>>> threshold=200<br>
>>>>> ><br>
>>>>><br>
>>>>> Here you chose a threshold of 200. That's 200 cells, so about 200 X<br>
>>>>> 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)<br>
>>>>> you<br>
>>>>> are getting over 19,000 tiny little catchments. Are you sure that's<br>
>>>>> what<br>
>>>>> you want??<br>
>>>>><br>
>>>>> Finally, you're trying to get raster values for 19,000 tiny vector<br>
>>>>> 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<br>
>>>>> causing<br>
>>>>> the NULL value error.<br>
>>>>><br>
>>>>> In short: I think you'll need to match the resolution of the DEM to<br>
>>>>> that<br>
>>>>> of the rainfall data. If the rainfall is only as accurate as 1 data<br>
>>>>> value<br>
>>>>> per 730 <a href="http://sq.km" target="_blank">sq.km</a>. then you will be able to do vector-raster analyses<br>
>>>>> only<br>
>>>>> at<br>
>>>>> that resolution = i.e. continent scale maps.<br>
>>>>><br>
>>>>> HTH<br>
>>>>> --<br>
>>>>> Micha<br>
>>>>><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<br>
>>>>> 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<br>
>>>>> '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<br>
>>>>> 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<br>
>>>>> 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<br>
>>>>> 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<br>
>>>>> 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<br>
>>>>> file,<br>
>>>>> >> so<br>
>>>>> >> you could say it's "wrong" if you want to keep the original<br>
>>>>> 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>
>>>>><br>
>>>><br>
>>>><br>
>>>> --<br>
>>>> Kind Regards<br>
>>>> TS Gumede<br>
>>>> CSIR, Meraka Institute<br>
>>>> 072 258 1650<br>
>>>><br>
>>>><br>
>>><br>
>>><br>
>>> --<br>
>>> Kind Regards<br>
>>> TS Gumede<br>
>>> CSIR, Meraka Institute<br>
>>> 072 258 1650<br>
>>><br>
>>><br>
>><br>
>><br>
>> --<br>
>> Kind Regards<br>
>> TS Gumede<br>
>> CSIR, Meraka Institute<br>
>> 072 258 1650<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>-- <br>Kind Regards<br>TS Gumede<br>CSIR, Meraka Institute<br>072 258 1650<br><br>