June 28, 2008
GDAL raster data tips and tricks
Reading GRASS data through GDAL/OGR support
Example 1: We write out a GRASS raster map to GeoTIFF -- this format
includes the coordinates within the file's metadata:
gdal_translate -of Gtiff /usr/local/share/grassdata/spearfish/PERMANENT/cellhd/soils soilmap.tif
ogr2ogr roadsmap.shp /usr/local/share/grassdata/spearfish/PERMANENT/vector/roads/head
Fast image display with tiling
If you want fast access you might want to try converting e.g. a BIL files to a tiled TIFF, and build overviews. You can build overviews for BIL too, but it can't be directly tiled:
# add -co "PROFILE=BASELINE" for TIF/TFW
gdal_translate source_bil global30.tif -co "TILED=YES" -co "TFW=YES" -co "PROFILE=BASELINE"
gdaladdo global30.tif 2 4 8 16
GDAL performance problem?
GDAL_CACHEMAX is normally a number of megabytes (default is 10 or so). So something like:
gdal_translate -of GTIFF -co TILED=YES --config GDAL_CACHEMAX 120 madison_1f_01.jpg madison_1f_01.tif
would use a 120MB cache.
GDAL and 1 bit maps
With a trick you can get those:
gdal_merge.py -co NBITS=1 -o dst.tif src.tif
Generate 8 bit maps for Mapserver
gdal_translate -scale in.tif out.tif
Note: As of MapServer 4.4 support has been added for classifying non-8bit raster inputs
A "proper" conversion would involve a colorspace transformation on the RGB image into IHS or something like that, and then taking the intensity. GRASS can do things like that.
Generate an OGC WKT (SRS)
In WKT the ellipsoid is described by two parameters: the semi-major axis and the inverse flattening. For a sphere the flattening is 0 and so the inverse flattening is infinity.
# in the GDAL source code:
./testepsg '+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m'
WKT[+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m] =
GEOGCS["GRS 1980(IUGG, 1980)",
Simplified WKT[+proj=lcc +lat_1=35 +lat_2=65 +lat_0=52 +lon_0=10 +x_0=4000000 +y_0=2800000 +ellps=GRS80 +units=m] =
GEOGCS["GRS 1980(IUGG, 1980)",
Extracting spatial subset (subregion)
W N E S
gdal_translate -of GTiff -projwin 636861 5152686 745617 5054047.5 \
Fixing broken projection/datum info for raster data
gdal_translate -of HFA -a_srs epsg:32735 /cdrom/173072lsat.img \
# or, using a WKT file
gdal_translate -of HFA -a_srs file.prj /cdrom/173072lsat.img \
Merge various import maps, re-project on the fly and extract spatial subset according to current GRASS region
eval `g.region -g`
gdalwarp -te $w $s $e $n *.TIF \
Export to (limited) TIFF readers such as ArcView, or ImageMagick
Many tools have trouble reading multi-band TIFFs with "band interleaving", the GDAL output default. Best is to use the INTERLEAVE=PIXEL creation option. Just add to the gdal_translate command line:
Inserting metadata (metadata tags)
gdal_translate -outsize 37.5% 37.5% \
-mo TIFFTAG_XRESOLUTION=300 -mo TIFFTAG_YRESOLUTION=300 \
Raster map reprojection (warping)
gdalwarp -t_srs '+init=epsg:26591 +towgs84=-225,-65,9' test1.tif \
Raster map reprojection (warping) maintaining NULL values (sea etc):
gdalwarp -r bilinear -tr 1000 1000 \
-srcnodata "-32768" -dstnodata "-32768" \
-wo "INIT_DEST=-32768" \
-t_srs epsg:32632 italy_LL.tif italy_UTM32.tif
Reprojecting external map to current GRASS location externally
gdalwarp -t_srs "`g.proj -wf`" aster.tif aster_tmerc.tif
Cut out region of interest with gdalwarp (in target coords)
Add to command line (insert values instead of letters of course:
#damn order, differs from -projwin!!
-te W S E N
Merging many small adjacent DEMs into one big map (A)
This needs GDAL compiled with Python and numpy installed:
# if not installed in standard site-packages directory
gdal_merge.py -v -o spearfishdem.tif -n "-32768" d*.tif
Merging many small adjacent DEMs into one big map (B)
Even easier, just use gdalwarp:
gdalwarp C_1mX1m/dtm*.tif big.tif
Or just a few tiles:
gdalwarp C_1mX1m/dtm0010[4-5]* big_selection.tif
Merge various map/bands into a RGB composite
gdal_merge.py -of HFA -separate band1.img band2.img band3.img -o out.img
GDAL gdalwarp interpolation comments
Which method -rn, rb, -rc or -rcs should one use for DEM and which for data like e.g. Landsat TM reprojecting?
-tps: Enable use of thin plate spline transformer based on available GCPs.
-rn: Use nearest neighbour resampling (default, fastest algorithm, worst interpolation quality).
-rb: Use bilinear resampling.
-rc: Use cubic resampling.
-rcs: Use cubic spline resampling (slowest algorithm).
I would suggest -rb for DEMs, and one of the cubic kernels for landsat data. Of course, there are various factors that you should take into account. Using -rb (bilinear) for the DEM will perform local averaging of the nearby pixel values in the source. This give reasonable results without introducing any risky "overshoot" effects you might see with cubic that could be disturbing for analysis or visualization in a DEM. The cubic should in theory do better at preserving edges and general visual crispness than using bilinar or nearest neighbour. However, if you are wanting to do analysis with the landsat (such as multispectral classification) I would suggest just using -rn (nearest neighbour) so as to avoid causing odd effects to the spectral values.
Nobody can't tell you what method should be used in your case. Generally speaking, in the case of upsampling spline and cubic interpolators are more suitable (-rcs and -rc). In the case of downsampling and the same resolution it is completely up to you what method looks better. Just try them all and select the one which is most appropriate for you.
Geocoding with 'gdal_translate'
As far as I know there is not on-screen method for doing this, but it certainly isn't too difficult with a little bit of semi-manual work. Open two OpenEV views, one with the unreferenced image, one with the geo-reference base you want to use. Move your cursor over the non-referenced one (let's call it image1), record (read: write down!) the pixel x,y values. Then look at the same location in image2. Write down the geocoordinate for the pixel. You should have four numbers for each location you want to pin the image to. And so on and so on. Then use gdal_translate to translate image1.tif to image1_georefd.tif but adding the -GCP parameter for each set of coordinates. Like so...
gdal_translate -gcp 1 1 500000 5000000 \
-gcp 200 400 550000 5250000 image1.tif \
Reading HDF ASTER
To select a channel and warp to UTM (or whatever is inside):
gdalwarp HDF4_SDS:ASTER_L1B:"pg-PR1B0000-2002031402_100_001":2 aster_2.tif