Используйте ccrs.epsg () для построения шейп-файла периметра почтового индекса в системе координат EPSG 4326

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

Согласно источнику, эти данные находятся в системе координат EPSG 4326. Когда я пытаюсь построить данные

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.io.img_tiles as cimgt 
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

# Create a Stamen terrain background instance
stamen_terrain = cimgt.Stamen('terrain-background')
fig = plt.figure(figsize = (mapsize,mapsize))
ax = fig.add_subplot(1, 1, 1, projection=stamen_terrain.crs)

# Set range of map, stipulate zoom level
ax.set_extent([-122.7, -121.5, 37.15, 38.15], crs=ccrs.Geodetic())
ax.add_image(stamen_terrain, 8, zorder = 0)

# Add city borders
shape_feature = ShapelyFeature(Reader(shapefile).geometries(), ccrs.epsg(4326), 
                linewidth = 2, facecolor = (1, 1, 1, 0), 
                edgecolor = (0.3, 0.3, 0.3, 1))
ax.add_feature(shape_feature, zorder = 1)

Я вижу следующую ошибку:

ValueError: EPSG code does not define a projection

Возможно связано - когда я смотрю на функцию ccrs.epsg(), она говорит, что этот код EPSG не поддерживается

help(ccrs.epsg)
Help on function epsg in module cartopy.crs:

epsg(code)
    Return the projection which corresponds to the given EPSG code.

    The EPSG code must correspond to a "projected coordinate system",
    so EPSG codes such as 4326 (WGS-84) which define a "geodetic coordinate
    system" will not work.

    Note
    ----
        The conversion is performed by querying https://epsg.io/ so a
        live internet connection is required.

Учитывая этот результат, я также попытался использовать ccrs.Geodetic():

# Add city borders
shape_feature = ShapelyFeature(Reader(shapefile).geometries(), ccrs.Geodetic(), 
                linewidth = 2, facecolor = (1, 1, 1, 0), 
                edgecolor = (0.3, 0.3, 0.3, 1))
ax.add_feature(shape_feature, zorder = 1)

Но при этом также не отображаются периметры почтового индекса и отображается это предупреждающее сообщение:

UserWarning: Approximating coordinate system <cartopy._crs.Geodetic object at 0x1a2d2375c8> with the PlateCarree projection.
  'PlateCarree projection.'.format(crs))

Я тоже пробовал ccrs.PlateCarree(), но безуспешно. Пожалуйста помоги!


python matplotlib shapefile cartopy shapely
person Michael Boles    schedule 14.07.2019    source источник
comment
Исходя из содержания ZIPCODE.prj, проекция равна NAD_1983_UTM_Zone_10N. Так что это не EPSG 4326. Подробнее см. epsg.io/26910.   -  person swatchai    schedule 15.07.2019
comment
Спасибо. Еще пробовал ccrs.epsg(26910) - ошибки не возвращает, но и границы не рисует - есть мысли, что дальше попробовать?   -  person Michael Boles    schedule 15.07.2019

Ответы (1)


arrow_upward
2
arrow_downward

Чтобы отобразить данные из разных источников вместе, нужно объявить правильный coordinate reference system для каждого набора данных. В случае с шейп-файлом вы можете найти его в прилагаемом к нему xxx.prj файле.

Вот рабочий код:

import cartopy.io.img_tiles as cimgt 
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

shapefile_name= "./data/ZIPCODE.shp"
mapwidth, mapheight = 8, 8
pad = 0.25

stamen_terrain = cimgt.Stamen('terrain-background')
stm_crs = stamen_terrain.crs

fig = plt.figure(figsize = (mapwidth, mapheight))
ax = fig.add_subplot(1, 1, 1, projection=stm_crs)  #world mercator

# Set extent of map
ax.set_extent([-123.3-pad, -121.5+pad, 37.05-pad, 38.75+pad], crs=ccrs.Geodetic())
# Plot base map
ax.add_image(stamen_terrain, 8, zorder=0)

# Add polygons from shapefile
# Note: the use of `ccrs.epsg(26910)`
shape_feature = ShapelyFeature(Reader(shapefile_name).geometries(), ccrs.epsg(26910))

# You can choose one of the 2 possible methods to plot
# ... the geometries from shapefile
# Styling is done here.
method = 1
if method==1:
    # iteration is hidden
    ax.add_feature(shape_feature, facecolor='b', edgecolor='red', alpha=0.4, zorder = 15)
if method==2:
    # iterate and use `.add_geometries()`
    # more flexible to manipulate particular items
    for geom in shape_feature.geometries():
        ax.add_geometries([geom], crs=shape_feature.crs, facecolor='b', edgecolor='red', alpha=0.4)

plt.show()

Выходной график:

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

person swatchai    schedule 18.07.2019
comment
Ваш method==2 пример весьма полезен. Благодаря этому я смог раскрасить полигоны по значению, заданному в отдельном фрейме данных: for counter, geom in enumerate(shape_feature.geometries()): if data['Area'][counter] < 50: if data['Population'][counter] > 500: ax.add_geometries([geom], crs=shape_feature.crs, facecolor=color[counter], edgecolor='k', alpha=0.8) else: continue ` - person Michael Boles; 27.07.2019
comment
@MichaelBoles :) - person swatchai; 28.07.2019