I have a bit of code from here that I am using for plotting a geotiff image (Landsat):
from mpl_toolkits.basemap import Basemap from osgeo import osr, gdal import matplotlib.pyplot as plt import numpy as np # Read the data and metadata ds = gdal.Open(filename) data = ds.ReadAsArray() gt = ds.GetGeoTransform() proj = ds.GetProjection() xres = gt[1] yres = gt[5] # get the edge coordinates and add half the resolution # to go to center coordinates xmin = gt[0] + xres * 0.5 xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5 ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5 ymax = gt[3] - yres * 0.5 ds = None # create a grid of xy coordinates in the original projection xy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres] # Create the figure and basemap object m = Basemap(projection='robin', lon_0=0, resolution='c') # Create the projection objects for the convertion inproj = osr.SpatialReference() inproj.ImportFromWkt(proj) # Get the target projection from the basemap object outproj = osr.SpatialReference() outproj.ImportFromProj4(m.proj4string) size = xy_source[0,:,:].size ct = osr.CoordinateTransformation(inproj, outproj) xy_target = np.array(ct.TransformPoints(xy_source.reshape(2, size).T))
But, it fails at ct.TransformPoints(xy_source.reshape(2, size).T))
and I’m not sure why. The error it gives me:
TypeError: in method ‘CoordinateTransformation_TransformPoints’, argument 1 of type ‘OSRCoordinateTransformationShadow *’
Which I do not understand. Any OSR guru’s out there?
Thanks for reading.
EDIT 1 The projection of my .TIFF
>>> print proj Out[20]: 'PROJCS["WGS 84 / UTM zone 34N",GEOGCS["WGS 84",DATUM["WGS_1984", SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]], AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0], UNIT["degree",0.0174532925199433],AUTHORITY["EPSG","4326"]], PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",21],PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000],PARAMETER["false_northing",0], UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32634"]]'
Also,
>>> m.proj4string Out[43]: '+lon_0=0.0 +y_0=8615499.05007 +R=6370997.0 +proj=robin +x_0=16986796.165 +units=m '
Advertisement
Answer
The solution is to install the proj4
package. Otherwise, the input projection is not understood by osr
.
conda install -c https://conda.binstar.org/jgomezdans proj4