如何使用matplotlib-basemap正確地投影tif圖像

Question

我嘗試使用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' 和 '不'， '之前ano_DOY ..' 'ano_DOY ..' 之間）

Answer 1

我我不確定你在做什麼錯你自己的插入/重新投影，但它可以做得更簡單。

contourf接受latlon關鍵字，如果爲true，則接受緯度/經度輸入並自動將其轉換爲地圖投影。所以：

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定義。