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 <[email protected]> > > > 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 <[email protected]> >> > 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 >> >> <[email protected]>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 >> >>> <[email protected]>wrote: >> >>> >> >>>> Hi >> >>>> >> >>>> Which module do I use to change the resolutions? >> >>>> >> >>>> >> >>>> 2010/7/6 <[email protected]> >> >>>> >> >>>>> 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 <[email protected]> >> >>>>> > >> >>>>> >> > 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 > > -- Kind Regards TS Gumede CSIR, Meraka Institute 072 258 1650
_______________________________________________ grass-user mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-user
