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
_______________________________________________ grass-user mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-user
