Как правильно проецировать изображение tif с помощью matplotlib-basemap

Я пытаюсь отобразить растровое изображение, используя gdal и matplotlib-basemap.

Здесь я объясняю свою попытку использовать функцию basemap.interp, для полного структурированного обзора моего процесса, пожалуйста, посмотрите мой Блокнот IPython. Сначала мой код для загрузки и проецирования растра.

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

Это построит следующее:

данные смещены по сравнению с береговой линией

Особенно хорошо видно, что на восточном побережье Северной и Южной Америки наблюдается смещение растровых данных и береговых линий.

Я понятия не имею, как адаптировать свой код, чтобы мои данные были преобразованы в правильную проекцию.

Для чего это стоит: Мой используемый растровый файл tif (если вы загружаете его, ставите «-» между «a» и «no», перед «ano_DOY..» после «a-no_DOY..»)


python matplotlib matplotlib-basemap gdal
person Mattijn    schedule 27.09.2013    source источник

Ответы (1)


arrow_upward
3
arrow_downward

Я не уверен, что вы делаете неправильно с вашей собственной интерполяцией/перепроецированием, но это можно сделать еще проще.

contourf принимает ключевое слово latlon, которое, если оно истинно, принимает входные данные широты и долготы и автоматически преобразует их в проекцию карты. Так:

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, чтобы удалить отражение вашего входного растра, просто личное предпочтение.

person Rutger Kassies    schedule 27.09.2013
comment
Спасибо за Ваш ответ. Это выглядит здорово! Но.. Я не могу воспроизвести ваши результаты.. Моя карта остается белой, за исключением части «добавить дополнительную информацию». Какую версию Basemap вы используете? Я установил версию 1.02 из дополнительных плагинов с сайта pythonxy. - person Mattijn; 27.09.2013
comment
Ключевое слово latlon введено в версии 1.05 (я использую 1.06). Я обновил свой ответ, чтобы выполнить «ручное» преобразование координат: xx, yy = m(lons,lats). Это работает для вас? - person Rutger Kassies; 27.09.2013
comment
Да! Это работает идеально. В настоящее время буду использовать это и скоро обновлю до более новой версии, так как 1.02 сильно устарела, я вижу - person Mattijn; 27.09.2013