> 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 <akasand...@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 >> <akasand...@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=...@permanent drainage=flow_direction >>> basin=catch accumulation=acc threshold=1 memory=300 stream=str >>> >>> # convert catch raster to polygon vector >>> r.to.vect in=ca...@permanent out=catchments feature=area >>> >>> g.region rast=rainfall -p >>> >>> # Calculate univariate statistics >>> v.rast.stats -c vector=catchme...@permanent >>> raster=rainf...@permanentcolpre=precp >>> >>> >>> >>> On Wed, Jul 7, 2010 at 9:18 AM, Sandile Gumede >>> <akasand...@gmail.com>wrote: >>> >>>> Hi >>>> >>>> Which module do I use to change the resolutions? >>>> >>>> >>>> 2010/7/6 <mi...@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 <mi...@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. > > > _______________________________________________ grass-user mailing list grass-user@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/grass-user