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=dem@PERMANENT drainage=flow_direction basin=catch accumulation=acc threshold=1 memory=300 stream=str
convert catch raster to polygon vector
r.to.vect in=catch@PERMANENT out=catchments feature=area
g.region rast=rainfall -p
Calculate univariate statistics
v.rast.stats -c vector=catchments@PERMANENT raster=rainfall@PERMANENT colpre=precp
On Wed, Jul 7, 2010 at 9:18 AM, Sandile Gumede <akasandile@gmail.com> wrote:
Hi
Which module do I use to change the resolutions?
2010/7/6 <micha@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=rainfallProjection of input dataset and current location appear to match
100%
r.in.gdal complete. Raster map 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: 3Here, 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=SRTMDEMProjection of input dataset and current location appear to match
100%
r.in.gdal complete. Raster map 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: 420894Your 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=200Here 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 …
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=precipDBMI-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.tifIs there something I have done wrong in these steps or there is something
wrong with my coordinates?2010/7/6 <micha@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 1650This mail was received via Mail-SeCure System.
–
Kind Regards
TS Gumede
CSIR, Meraka Institute
072 258 1650This 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