График GDAL с использованием matplotlib Базовая карта

Я хотел бы построить растровый tiff ( скачать -723Kb) с использованием matplotlib Basemap. Координаты проекции растра находятся в метрах:

In [2]: path = r'albers_5km.tif' raster = gdal.Open(path, gdal.GA_ReadOnly) array = raster.GetRasterBand(20).ReadAsArray() print ('Raster Projection:\n', raster.GetProjection()) print ('Raster GeoTransform:\n', raster.GetGeoTransform()) Out [2]: Raster Projection: PROJCS["unnamed",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["Albers_Conic_Equal_Area"],PARAMETER["standard_parallel_1",15],PARAMETER["standard_parallel_2",65],PARAMETER["latitude_of_center",30],PARAMETER["longitude_of_center",95],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]]] Raster GeoTransform: (190425.8243, 5000.0, 0.0, 1500257.0112, 0.0, -5000.0) 

Если я попытаюсь построить это с использованием проекции Робин, используя contourf с latlon=False чем x и y, предполагается, что это координаты проекции карты (см. Документы , я думаю, это то, что у меня есть).

Но если я посмотрю на сюжет, я заметил, что он расположен внизу слева очень маленьким:

Нижняя левая

Используя этот код:

 In [3]: xy = raster.GetGeoTransform() x = raster.RasterXSize y = raster.RasterYSize lon_start = xy[0] lon_stop = x*xy[1]+xy[0] lon_step = xy[1] lat_start = xy[3] lat_stop = y*xy[5]+xy[3] lat_step = xy[5] fig = plt.figure(figsize=(16,10)) map = Basemap(projection='robin',resolution='c',lat_0=0,lon_0=0) lons = np.arange(lon_start, lon_stop, lon_step) lats = np.arange(lat_start, lat_stop, lat_step) xx, yy = np.meshgrid(lons,lats) levels = [array.min(),-0.128305,array.max()] map.contourf(xx, yy,array, levels, cmap=cm.RdBu_r, latlon=False) map.colorbar(cntr,location='right',pad='10%') map.drawcoastlines(linewidth=.5) map.drawcountries(color='red') 

В конце концов я не хочу иметь мировоззрение, но подробный обзор. Но это дает мне уровень масштабирования, когда береговые линии и страны рисуются, но данные снова помещаются в нижний левый угол, но не так малы, как в предыдущее время:

снова внизу слева

Используя следующий код:

 In [4]: extent = [ xy[0],xy[0]+x*xy[1], xy[3],xy[3]+y*xy[5]] width_x = (extent[1]-extent[0])*10 height_y = (extent[2]-extent[3])*10 fig = plt.figure(figsize=(16,10)) map = Basemap(projection='stere', resolution='c', width = width_x , height = height_y, lat_0=40.2,lon_0=99.6,) xx, yy = np.meshgrid(lons,lats) levels = [array.min(),-0.128305,array.max()] map.contourf(xx, yy, array, levels, cmap=cm.RdBu_r, latlon=False) map.drawcoastlines(linewidth=.5) map.drawcountries(color='red') 

One Solution collect form web for “График GDAL с использованием matplotlib Базовая карта”

Вы можете использовать следующий код для преобразования координат, он автоматически принимает проекцию из вашего растра в качестве источника и проекции из вашего объекта Basemap в качестве целевой системы координат.

импорт

 from mpl_toolkits.basemap import Basemap import osr, gdal import matplotlib.pyplot as plt import numpy as np 

Преобразование координат

 def convertXY(xy_source, inproj, outproj): # function to convert coordinates shape = xy_source[0,:,:].shape size = xy_source[0,:,:].size # the ct object takes and returns pairs of x,y, not 2d grids # so the the grid needs to be reshaped (flattened) and back. ct = osr.CoordinateTransformation(inproj, outproj) xy_target = np.array(ct.TransformPoints(xy_source.reshape(2, size).T)) xx = xy_target[:,0].reshape(shape) yy = xy_target[:,1].reshape(shape) return xx, yy 

Чтение и обработка данных

 # Read the data and metadata ds = gdal.Open(r'albers_5km.tif') 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 fig = plt.figure(figsize=(12, 6)) m = Basemap(projection='robin', lon_0=0, resolution='c') # Create the projection objects for the convertion # original (Albers) inproj = osr.SpatialReference() inproj.ImportFromWkt(proj) # Get the target projection from the basemap object outproj = osr.SpatialReference() outproj.ImportFromProj4(m.proj4string) # Convert from source projection to basemap projection xx, yy = convertXY(xy_source, inproj, outproj) # plot the data (first layer) im1 = m.pcolormesh(xx, yy, data[0,:,:].T, cmap=plt.cm.jet) # annotate m.drawcountries() m.drawcoastlines(linewidth=.5) plt.savefig('world.png',dpi=75) 

введите описание изображения здесь

Если вам нужно, чтобы пиксельное местоположение было на 100% правильным, вы можете проверить создание массивов координат, действительно заботясь о себе (потому что я вообще не был). Этот пример, надеюсь, поставит вас на правильный путь.

  • IndexError с Basemap.contour () при использовании определенных прогнозов
  • наложение базовой карты на контуры
  • Базовая карта с Python 3.5 Anaconda в Windows
  • Python колчан и pcolormesh не выстраиваются точно вправо
  • Как показать метку шейп-файла в легенде базовой карты python?
  • Python: скопируйте базовую карту или удалите данные с рисунка
  • Построение океанов на картах с использованием базовой карты и питона
  • Базовая карта Matplotlib: всплывающее окно
  •  
    Interesting Posts for Van-Lav

    Преобразование объекта строки sqlalchemy в python dict

    Переопределить метод на уровне экземпляра

    Начать скрининг по маршруту Flask

    Объединение файлов hdf5

    Как преобразовать мой базовый код TensorFlow на основе подачи, чтобы использовать «Набор данных»?

    Префикс prepend для перечисления элементов со списком

    gevent: недостаток в нерестах большого количества зеленых?

    как узнать, удаляется ли объект в python

    Как использовать сложные SQL-скрипты с python

    django querysets + memcached: лучшие практики

    Как я могу обойти ошибку модуля kivy: ImportError: Ошибка загрузки DLL: указанного модуля не удалось найти?

    Форматирование строки с форматом «{0: d}» дает код неизвестного формата «d» для объекта типа «float»

    TypeError: объект 'int' не повторяется в python

    Как организовать три списка таким образом, чтобы сумма соответствующих элементов, если они больше, появляются сначала?

    Изменение моего пути python: helloworld.py возвращает команду не найденной –

    Python - лучший язык программирования в мире.