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') 

이

은 다음과 플롯됩니다

이

북미와 남미의 동부 해안에 래스터 데이터와 해안선의 오프셋 (offset)가있는 것을 볼 특히 분명하다.

내 데이터를 올바른 방법으로 변형되도록 코드를 수정하는 방법이 단서가 아닙니다. 그것은 가치가 무엇인지에 대한

: My used raster tif file (당신은 그것을두고 다운로드 할 경우 '-'이전에 'A'와 '아니오', 'ano_DOY ..'후 'ano_DOY ..'사이) 나 '

Answer 1

자신의 내삽/재 투영으로 무엇을 잘못하고 있는지 확신 할 수 없지만 더 간단하게 수행 할 수 있습니다.

contourf은 latlon 키워드를 허용하며, 사실 경우 위도/경도 입력을 받아 들여지도 투영으로 자동 변환합니다. 그래서 : 나는 당신의 입력 래스터의 반전, 단지 개인적인 취향을 제거하기 위해 lats 정의를 변경

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') 

참고.