20090327

Georeferencing raster data

Two maps:
  • "Soils of Arusha and Monduli Area 1:250.000" (Soil Research Institute, Ottawa, Canada, 1974)
  • "Tanzania Eco-Climatic Zones" (Institute of Resource Assestment).
were georeferenced in WGS84LL using Georeferencing Qgis Plugin.

Other sources can be georeferenced using gdal utilities.
For testing purposes is possible to download high resolution imagery covering part of Area of Intervention from Google dataset using GoogleMV. Level 17 has enough resolution.
The result is a big *.jpg file and a *.map file (Ozi Explorer format, see here info on the format).
The *.map is simple text, and coordinates of the corners of the image in latlongwgs84 can be found near the end of the file:
MMPLL,1, 36.754167, -3.076833
MMPLL,2, 36.903752, -3.076833
MMPLL,3, 36.903752, -3.233316
MMPLL,4, 36.754167, -3.233316
(in red are marked the Upper Left and Lower Right corners coordinates).

gdal_translate can be used now:

gdal_translate -a_srs EPSG:4326 -a_ullr 36.754167 -3.076833 36.903752 -3.233316 -co "COMPRESS=JPEG" -co "PHOTOMETRIC=YCBCR" -co "JPEG_QUALITY=100" source.jpg target.tif

  • the -a_srs parameter assign an EPSG code to the image
  • the -a_ullr parameter define upper left and lower right corners coordinates (taken from the *.map file)
  • the -co parameters define the compression algoritm used
The result is a GeoTIFF with this header (read it with gdalinfo):

Driver: GTiff/GeoTIFF
Files: target.tif
Size is 13943, 14608
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.2572235630016,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (36.754167000000002,-3.076833000000000)
Pixel Size = (0.000010728322456,-0.000010712144031)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  SOURCE_COLOR_SPACE=YCbCr
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  36.7541670,  -3.0768330) ( 36d45'15.00"E,  3d 4'36.60"S)
Lower Left  (  36.7541670,  -3.2333160) ( 36d45'15.00"E,  3d13'59.94"S)
Upper Right (  36.9037520,  -3.0768330) ( 36d54'13.51"E,  3d 4'36.60"S)
Lower Right (  36.9037520,  -3.2333160) ( 36d54'13.51"E,  3d13'59.94"S)
Center      (  36.8289595,  -3.1550745) ( 36d49'44.25"E,  3d 9'18.27"S)
Band 1 Block=13943x16 Type=Byte, ColorInterp=Red
Band 2 Block=13943x16 Type=Byte, ColorInterp=Green
Band 3 Block=13943x16 Type=Byte, ColorInterp=Blue

The file is huge, so pyramid layers can be added to fasten the opening in GIS with gdaladdo:

gdaladdo --config COMPRESS_OVERVIEW JPEG --config PHOTOMETRIC_OVERVIEW YCBCR --config INTERLEAVE_OVERVIEW PIXEL target.tif 2 4 8 16

In this way the last lines of the header of GeoTIFF are substituted with:

Band 1 Block=13943x16 Type=Byte, ColorInterp=Red
  Overviews: 6972x7304, 3486x3652, 1743x1826, 872x913
Band 2 Block=13943x16 Type=Byte, ColorInterp=Green
  Overviews: 6972x7304, 3486x3652, 1743x1826, 872x913
Band 3 Block=13943x16 Type=Byte, ColorInterp=Blue
  Overviews: 6972x7304, 3486x3652, 1743x1826, 872x913
You can obtain a shapefile with the coverage of the GeoTiff with:
  • gdaltindex target.shp target.tif
Links:

GoogleMV:
GDAL Utilities:
OZI *.map format:




Powered by ScribeFire.

No comments:

Post a Comment