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