2

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.

2 Answers 2

1

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.

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

Try to go to spatial reference

To find specific proj4 definition

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.