Reprojecting data

Another map has been were georeferenced in WGS84LL using Georeferencing Qgis Plugin:

• "Agroecological zones of Tanzania", Ministry Of Agriculture, Food and Cooperative

Some reprojection has been done on raster maps using gdal tools:

No Datum change:

gdalwarp -t_srs EPSG:32737 4326_HRImage.tif 32737_HRImage.tif

move the source WGS84LL source GeoTiff file to the WGS84UTM37 huge, uncompressed target file, so you need to reduce it with

gdal_translate -co "COMPRESS=JPEG" -co "PHOTOMETRIC=YCBCR" -co "JPEG_QUALITY=100" 32737_HRImage.tif 32737_HRImage2.tif

and add pyramid overviews (load very fast in Qgis) with

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

Datum change:

gdalwarp -t_srs EPSG:32737 -s_srs "1285_TZ.prf" 21037_Topo.tif 32737_Topo.tif

move the source ARC60UTM37 source GeoTiff file to the WGS84UTM37 file, using the right transformation parameters written in:

1285_TZ.prf
+proj=utm +zone=37 +south +a=6378249.145 +b=6356514.96582849 +units=m +towgs84=-175,-23,-303

The process works backward: from standard EPSG_code --> to custom_CRS

In this case you can add right info to the header of the result with:

gdal_translate -a_srs "21037_WK.txt" source.tif 21037_target.tif

with 21037_WK.txt:

PROJCS[" Projection Name = Arc_1960_UTM_Zone_37S Units = meters GeoTIFF Units = meters",
    GEOGCS["Arc 1960",
        DATUM["Arc_1960",
            SPHEROID["Clarke 1880 (RGS)",6378249.145,293.4649999999983,
                AUTHORITY["EPSG","7012"]],
            AUTHORITY["EPSG","6210"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4210"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",39],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",10000000],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]

that can be extract with gdalinfo>info.txt on the original toposheets (see Defining the Area of Intervention)

In this way we obtained:

./Raster/21037:
21037_Agro_Eco_Zones.tif
21037_Arusha_soil.tif
21037_EcoClimatic_Zones.tif
21037_HRImage.tif
21037_Topo.tif

The same structure is repeated for 32737 and 4326.

With

gdaltindex index.shp *.tif

we obtained the coverage of all the rasters in every folder.

The same transformation procedures work with vector data:

ogr2ogr 4326_Index.shp 21037_Index.shp -t_srs EPSG:4326 -s_srs "1285.prf"

or

ogr2ogr mosaic2.shp mosaic.shp -t_srs EPSG:4326 -s_srs "+init=EPSG:21037 +wgs84=-175,-23,-303"

Links:
http://www.gdal.org/gdalwarp.html

Powered by ScribeFire.