Построение замаскированной Антарктиды с помощью шейп-файла или геопанды - PullRequest
0 голосов
/ 16 мая 2018

Я пытаюсь нанести данные вокруг Антарктиды, маскируя континент.Хотя я использую basemap и у него есть возможность легко замаскировать континенты, используя map.fillcontinents(), континент, рассматриваемый basemap, включает в себя шельфовые ледники, которые я не хочу маскировать.

Я пытался использовать geopandas из кода, найденного в Интернете.Это работает, за исключением того, что береговая линия создает нежелательную линию, которая, как я полагаю, является началом / концом многоугольника для Антарктиды:

import numpy as np
from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
from matplotlib.collections import PatchCollection
import geopandas as gpd
import shapely
from descartes import PolygonPatch

lats = np.arange(-90,-59,1)
lons = np.arange(0,361,1)
X, Y = np.meshgrid(lons, lats)
data = np.random.rand(len(lats),len(lons))

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))

fig=plt.figure(dpi=150)
ax = fig.add_subplot(111)

m = Basemap(projection='spstere',boundinglat=-60,lon_0=180,resolution='i',round=True)  
xi, yi = m(X,Y)

cf = m.contourf(xi,yi,data)

patches = []
selection = world[world.name == 'Antarctica']
for poly in selection.geometry:
    if poly.geom_type == 'Polygon':
        mpoly = shapely.ops.transform(m, poly)
        patches.append(PolygonPatch(mpoly))
    elif poly.geom_type == 'MultiPolygon':
        for subpoly in poly:
            mpoly = shapely.ops.transform(m, poly)
            patches.append(PolygonPatch(mpoly))
    else:
        print(poly, 'blah')
ax.add_collection(PatchCollection(patches, match_original=True,color='w',edgecolor='k'))

Geopandas output

Эта же строка появляется, когда я пытаюсь использовать другие шейп-файлы, такие как land , который можно бесплатно загрузить с Natural Earth Data .Поэтому я отредактировал этот шейп-файл в QGIS, чтобы удалить границы Антарктиды.Теперь проблема в том, что я не знаю, как замаскировать все, что внутри шейп-файла (и не мог найти, как это сделать).Я также попытался объединить предыдущий код с geopandas, установив linewidth=0 и добавив сверху созданный мной шейп-файл.Проблема в том, что они не совсем одинаковы:

Bad borders when using geopandas and shapefile

Любое предложение о том, как замаскировать с помощью шейп-файла или с геопандами, но без линии?

Редактировать: Используя предыдущий ответ Томаса Кхуна с моим отредактированным шейп-файлом, можно получить хорошо замаскированные Антарктиду / континенты, но береговая линия выходит за круглые края карты:

outside edges

Я загрузил здесь отредактированный шейп-файл, который я использовал, но это Естественный файл данных Земли на 50 м * шейп-файл без линии.

1 Ответ

0 голосов
/ 19 мая 2018

Вот пример того, как добиться того, чего вы хотите. Я в основном следовал примеру Basemap , как работать с shapefiles, и добавил немного shapely magic , чтобы ограничить контуры границ карты. Обратите внимание, что я сначала попытался извлечь контур карты из ax.patches, но это почему-то не сработало, поэтому я определил круг с радиусом boundinglat и преобразовал его, используя функцию преобразования координат Базовой карты.

import numpy as np

import matplotlib.pyplot as plt
from mpl_toolkits.basemap import Basemap
from matplotlib.collections import PatchCollection
from matplotlib.patches import Polygon

import shapely
from shapely.geometry import Polygon as sPolygon

boundinglat = -40 
lats = np.arange(-90,boundinglat+1,1)
lons = np.arange(0,361,1)
X, Y = np.meshgrid(lons, lats)
data = np.random.rand(len(lats),len(lons))

fig, ax = plt.subplots(nrows=1, ncols=1, dpi=150)

m = Basemap(
    ax = ax,
    projection='spstere',boundinglat=boundinglat,lon_0=180,
    resolution='i',round=True
)  
xi, yi = m(X,Y)

cf = m.contourf(xi,yi,data)

#adjust the path to the shapefile here:
result = m.readshapefile(
    'shapefiles/AntarcticaWGS84_contorno', 'antarctica',
    zorder = 10, color = 'k', drawbounds = False)

#defining the outline of the map as shapely Polygon:
rim = [np.linspace(0,360,100),np.ones(100)*boundinglat,]
outline = sPolygon(np.asarray(m(rim[0],rim[1])).T)

#following Basemap tutorial for shapefiles
patches = []
for info, shape in zip(m.antarctica_info, m.antarctica):
    #instead of a matplotlib Polygon, create first a shapely Polygon
    poly = sPolygon(shape)
    #check if the Polygon, or parts of it are inside the map:
    if poly.intersects(outline):
        #if yes, cut and insert
        intersect = poly.intersection(outline)        
        verts = np.array(intersect.exterior.coords.xy)
        patches.append(Polygon(verts.T, True))

ax.add_collection(PatchCollection(
    patches, facecolor= 'w', edgecolor='k', linewidths=1., zorder=2
))

plt.show()

Результат выглядит так:

the result of the above code

Надеюсь, это поможет.

...