График 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% правильным, вы можете проверить создание массивов координат, действительно заботясь о себе (потому что я вообще не был). Этот пример, надеюсь, поставит вас на правильный путь.

  • Базовая карта Python «Утверждение не выполнено»
  • Можно ли контролировать ориентацию маркера matplotlib?
  • Найти более простой способ группировать 2-мерные данные рассеяния в данные массива сетки
  • Массовая анимация Matplotlib
  • Исключить белые грани в Matplotlib / Basemap pcolor plot
  • Базовая карта на поверхности куба matplotlib.mplot3d
  • как маскировать данные конкретного массива на основе шейп-файла
  • Резервная сетка от центральных координат к внешним (то есть угловым) координатам
  • Сохранение наложения карт между графиками в matplotlib
  • Построение кривой кривой в базе данных Python
  • Построение океанов на картах с использованием базовой карты и питона
  •  
    Interesting Posts for Van-Lav

    Найти все таблицы в html с помощью BeautifulSoup

    Каталог поиска для определенной строки

    «Правильный» способ добавления сценариев python в приложение, отличное от python

    Numpy уникальная 2D-подматрица

    Как интегрировать любой питон Python с API статуса фиксации GitHub?

    Как проверить, является ли объект итератором в Python?

    Как хранятся и вычисляются истории контроля версий?

    Передача аргумента в соединение PyQt SIGNAL

    Количество глаголов, существительных и других частей речи с помощью NLTK python

    Создание фильтра нижних частот в SciPy – методы и единицы понимания

    Почему Pylint дает ошибку E0702, поднимая NoneType, в этом выражении raise?

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

    Python в Windows – как ждать нескольких дочерних процессов?

    Python + SSH Password auth (нет внешних библиотек или общедоступных / закрытых ключей)?

    Шифрование Python с помощью PyCrypto AES

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