python之如何使用 matplotlib-basemap 正确投影 tif 图像

我尝试使用 gdal 和 matplotlib-basemap 显示光栅图像。

我在这里解释了我使用 basemap.interp 函数的尝试，有关我的过程的完整结构化概述，请查看我的 IPython Notebook . 首先是加载和投影栅格的代码。

# Load Raster 
pathToRaster = r'I:\Data\anomaly//ano_DOY2002170.tif' 
raster = gdal.Open(pathToRaster, gdal.GA_ReadOnly) 
array = raster.GetRasterBand(1).ReadAsArray() 
msk_array = np.ma.masked_equal(array, value = 65535) 
print 'Raster Projection:\n', raster.GetProjection() 
print 'Raster GeoTransform:\n', raster.GetGeoTransform() 
 
# Project raster image using Basemap and the basemap.interp function 
map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0) 
 
datain = np.flipud( msk_array ) 
 
nx = raster.RasterXSize 
ny = raster.RasterYSize 
 
xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on the grid 
yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on the grid 
 
lons = np.arange(-180,180,0.25) #from raster.GetGeoTransform() 
lats  = np.arange(-90,90,0.25)  
 
lons, lats = np.meshgrid(lons,lats)  
xout,yout = map(lons, lats) 
dataout = mpl_toolkits.basemap.interp(datain, xin, yin, xout, yout, order=1) 
 
levels = [-1000,-800,-600,-400,-200,0,200,400,600,800,1000] 
cntr = map.contourf(xout,yout,dataout, levels,cmap=cm.RdBu) 
cbar = map.colorbar(cntr,location='bottom',pad='15%') 
 
# Add some more info to the map 
cstl = map.drawcoastlines(linewidth=.5) 
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' )  
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey') 
boun = map.drawmapboundary(linewidth=0.5, color='grey') 

这将绘制以下内容:

特别明显的是，在北美和南美的东海岸，栅格数据与海岸线存在偏移。

我不知道如何调整我的代码，以便我的数据在正确的投影中进行转换。

物有所值:My used raster tif file (如果你下载它，它会在 'a' 和 'no' 之间放置一个 '-'，在 'ano_DOY..' 之前，在 'a-no_DOY..' 之后)

请您参考如下方法:

我不确定你自己的插值/重新投影做错了什么，但它可以更简单地完成。

contourf 接受一个 latlon 关键字，当该关键字为真时，接受纬度/经度输入并自动将其转换为 map 投影。所以:

datain = msk_array 
 
fig = plt.figure(figsize=(12,5)) 
map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0) 
 
ny, nx = datain.shape 
 
xin = np.linspace(map.xmin,map.xmax,nx) # nx is the number of x points on the grid 
yin = np.linspace(map.ymin,map.ymax,ny) # ny in the number of y points on the grid 
 
lons = np.arange(-180,180,0.25) #from raster.GetGeoTransform() 
lats  = np.arange(90,-90,-0.25)  
 
lons, lats = np.meshgrid(lons,lats) 
 
xx, yy = m(lons,lats) 
 
levels = [-1000,-800,-600,-400,-200,0,200,400,600,800,1000] 
cntr = map.contourf(xx, yy,datain, levels,cmap=cm.RdBu) 
 
cbar = map.colorbar(cntr,location='bottom',pad='15%') 
 
# Add some more info to the map 
cstl = map.drawcoastlines(linewidth=.5) 
meri = map.drawmeridians(np.arange(0,360,60), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey' )  
para = map.drawparallels(np.arange(-90,90,30), linewidth=.2, labels=[1,0,0,1], labelstyle='+/-', color='grey') 
boun = map.drawmapboundary(linewidth=0.5, color='grey') 

请注意，我更改了 lats 定义以删除输入栅格的翻转，这只是个人喜好。