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.