使用 matplotlib 底图绘制 .tif GDAL 栅格答案

【问题标题】：Plot .tif GDAL raster using matplotlib Basemap使用 matplotlib 底图绘制 .tif GDAL 栅格
【发布时间】：2014-10-11 14:08:10
【问题描述】：

我有一些来自 here 的代码用于绘制 geotiff 图像 (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))

但是，它在ct.TransformPoints(xy_source.reshape(2, size).T)) 失败了，我不知道为什么。它给我的错误：

TypeError：在“CoordinateTransformation_TransformPoints”方法中，“OSRCoordinateTransformationShadow *”类型的参数 1

我不明白。有没有 OSR 大师？

感谢阅读。

EDIT 1我的.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"]]'

还有，

>>> m.proj4string
Out[43]: '+lon_0=0.0 +y_0=8615499.05007 +R=6370997.0 +proj=robin +x_0=16986796.165 +units=m '

【问题讨论】：

听起来您的输入数据可能没有定义投影。 proj 是什么样的？（即在完成proj = ds.GetProjection() 之后输入print proj）。
尝试在顶部使用osr.UseExceptions()。 proj 和 m.proj4string 输入是什么？

标签： python linux matplotlib gdal matplotlib-basemap

【解决方案1】：

解决方案是安装proj4 包。否则，osr 无法理解输入投影。

conda install -c https://conda.binstar.org/jgomezdans proj4

【讨论】：