Is there any one who can help me? 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
_______________________________________________ grass-user mailing list [email protected] http://lists.osgeo.org/mailman/listinfo/grass-user
