[GRASS-user] Re: gdal_translate

Tue Jul 6 11:24:06 EDT 2010

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 <micha at 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.
>
>
>