How to define PROJ.4 string for GDAL projection

Question

I'm doing some zonal computations in GDAL using the code posted in the solution here: Issue Trying to create Zonal Statistics using Gdal and Python.

It runs fine for me using the example data posted with the question, but when I try with my own GeoTiff and shapefile, this step:

# Read projection info from input GeoTiff
raster_srs.ImportFromWkt(raster.GetProjectionRef())
target_ds.SetProjection(raster_srs.ExportToWkt())

# Rasterize zone polygon to raster
gdal.RasterizeLayer(target_ds, [1], lyr, burn_values=[1])

throws the following error:

ERROR 1: No PROJ.4 translation for destination SRS, coordinate
transformation initialization has failed.

gdal info:

Driver: GTiff/GeoTIFF
Files: data/ndvi/ndvi_c_dk_20150405.tif
Size is 29424, 30470
Coordinate System is:
LOCAL_CS["unnamed",
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (-19.000000000000000,37.999585000000003)
Pixel Size = (0.002413000000000,-0.002413000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  ( -19.0000000,  37.9995850) 
Lower Left  ( -19.0000000, -35.5245250) 
Upper Right (  52.0001120,  37.9995850) 
Lower Right (  52.0001120, -35.5245250) 
Center      (  16.5000560,   1.2375300) 
Band 1 Block=29424x1 Type=Byte, ColorInterp=Palette
  Color Table (RGB with 256 entries)

Knowing that the coordinate systems info is missing, I think I should explicitly define the Proj4 string. I've found examples of doing such like this:

wgs84_wkt = """
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.01745329251994328,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]"""

then I presume I could do

raster_srs.ImportFromWkt(wgs84_wkt)
target_ds.SetProjection(raster_srs.ExportToWkt())

But I have no idea how to define the correct Proj4 string in the first place.

score 1 · Accepted Answer · answered May 22 '15 at 07:01

If the dataset should cover the whole African continent, the CRS information is wrong, and WGS84 should be right.

If it is really a local CRS, covering 71m eastwards and 73.5m northwards, you have to set up a local CRS on the point of (0;0), which you have to collect with GPS or offical surveying information:

+proj=tmerc +lat_0=.... +lon_0=.... +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

For such a small area, projection distortions are very low, so you should get lucky with tmerc.

It is the whole African continent - thanks for explaining what the details were. I'm new to GDAL and GIS development. — user3591836, May 22 '15 at 15:06

score -1 · Answer 2 · answered May 21 '15 at 23:07

-1

Try to go to spatial reference

To find specific proj4 definition

answered May 21 '15 at 23:07

Patrick Sky

1

How to define PROJ.4 string for GDAL projection

2 Answers2