python - Plot .tif GDAL raster using matplotlib Basemap -
i have bit of code here using plotting geotiff image (landsat):
from mpl_toolkits.basemap import basemap osgeo import osr, gdal import matplotlib.pyplot plt import numpy np # read data , metadata ds = gdal.open(filename) data = ds.readasarray() gt = ds.getgeotransform() proj = ds.getprojection() xres = gt[1] yres = gt[5] # edge coordinates , add half resolution # go 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 grid of xy coordinates in original projection xy_source = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres] # create figure , basemap object m = basemap(projection='robin', lon_0=0, resolution='c') # create projection objects convertion inproj = osr.spatialreference() inproj.importfromwkt(proj) # target projection 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, fails @ ct.transformpoints(xy_source.reshape(2, size).t))
, i'm not sure why. error gives me:
typeerror: in method 'coordinatetransformation_transformpoints', argument 1 of type 'osrcoordinatetransformationshadow *'
which not understand. osr guru's out there?
thanks reading.
edit 1 projection of .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 '
the solution install proj4
package. otherwise, input projection not understood osr
.
conda install -c https://conda.binstar.org/jgomezdans proj4
Comments
Post a Comment