Совместите проекцию шейп-файла в картопии

Я пытаюсь создать карту Choropleth с использованием matplotlib и cartopy, для которой мне, очевидно, сначала нужно построить шейп-файл. Однако мне не удалось это сделать, хотя аналогичный вопрос был задан здесь и здесь. Я подозреваю, что неверно указаны либо проекция, либо границы.

В моем шейп-файле есть проекция

PROJCS["WGS_1984_UTM_Zone_32Nz",
    GEOGCS["GCS_WGS_1984",
        DATUM["WGS_1984",
            SPHEROID["WGS_84",6378137,298.257223563]],
        PRIMEM["Greenwich",0],
        UNIT["Degree",0.017453292519943295]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["False_Easting",32500000],
    PARAMETER["False_Northing",0],
    PARAMETER["Central_Meridian",9],
    PARAMETER["Scale_Factor",0.9996],
    PARAMETER["Latitude_Of_Origin",0],
    UNIT["Meter",1]]

и его можно загрузить здесь, где я говорю о vg250_2010-01-01.utm32w.shape.ebenen/vg250_ebenen-historisch/de1001/vg250_gem.shp

Мой код

#!/usr/local/bin/python
# -*- coding: utf8 -*-

import cartopy.crs as ccrs
import cartopy.io.shapereader as shpreader
import matplotlib.pyplot as plt

fname = 'path/vg250_gem.shp'
proj = ccrs.TransverseMercator(central_longitude=0.0,central_latitude=0.0,
                           false_easting=32500000.0,false_northing=0.0,
                           scale_factor=0.9996)
municipalities = list(shpreader.Reader(fname).geometries())
ax = plt.axes(projection=proj)
plt.title('Deutschland')
ax.add_geometries(municipalities,proj,edgecolor='black',facecolor='gray',alpha=0.5)
ax.set_extent([32458044.649189778*0.9, 5556418.748046352*1.1, 32465287.307457082*0.9, 5564153.5456742775*1.1],proj)
plt.show()

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

Traceback (most recent call last):
  File "***/src/analysis/test.py", line 16, in <module>
ax.set_extent([32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775],proj)
  File "/usr/local/lib/python2.7/site-packages/cartopy/mpl/geoaxes.py", line 652, in set_extent
ylim=self.projection.y_limits))
ValueError: Failed to determine the required bounds in projection coordinates. Check that the values provided are within the valid range (x_limits=[-20000000.0, 20000000.0], y_limits=[-10000000.0, 10000000.0]).
[Finished in 53.9s with exit code 1]

Для меня это не имеет смысла. Кроме того, экспериментируя с ccrs.UTM (), я получаю график, показывающий белую область. Буду признателен, если кто-нибудь скажет мне, как это исправить. Спасибо!


matplotlib cartopy projection
person Jhonny    schedule 08.10.2016    source источник
comment
Я ничего не знаю о картопии и т. Д., Но, читая предпоследнюю строку в сообщении об ошибке, есть некоторые обязательные ограничения по x и y, но ограничения в вашем вызове set_extent (вызов, вызвавший ошибку ) выходят за эти требуемые пределы.   -  person tsj    schedule 08.10.2016
comment
Это то, о чем я говорил. Это не имеет для меня смысла, поскольку я извлек границы шейп-файла с помощью fiona. Выше я пробовал поэкспериментировать с методом set_extent, в том числе полностью удалить эту строку, но это тоже не сработало.   -  person Jhonny    schedule 09.10.2016
comment
Независимо от того, какие границы вам были даны, они выходят за рамки, которые картография готова принять. Опять же, я не знаком с этими пакетами, но я заметил, что при построении proj вы предоставляете false_easting 32,5M. Предположительно это приведет к сдвигу в направлении x. В вашем обращении к set_extent у вас есть два значения около 32,5M, которые будут исправлены таким сдвигом. Глядя на документы для картографии, входные данные для set_extent - (x0, x1, y0, y1), но похоже, что вы могли указать (x0, y0, x1, y1), поэтому false_easting не вернет вас к приемлемому границы.   -  person tsj    schedule 10.10.2016

Ответы (1)


arrow_upward
2
arrow_downward

Я обнаружил две проблемы. Один из них - неправильное указание ограничений в вашем вызове set_extent, документация указывает, что [x0,x1,y0,y1] должен быть вводом, вы, кажется, дали [x0,y0,x1,y1]. Другая проблема, по-видимому, связана с ограничением картографии, насколько я могу судить. Похоже, что прогнозы, выходящие за пределы, указанные в сообщении об ошибке, всегда не удастся, и эти ограничения жестко запрограммированы. Вы можете просто отредактировать источник (эту строку в их последний выпуск), изменив -2e7 на -4e7, а также верхнюю границу. После этих исправлений ваш график создается без проблем: Deutschland Новая set_extent строка:

ax.set_extent([32458044.649189778*0.975, 32465287.307457082*1.025,5556418.748046352*0.9, 556415,3.5456742775*1.1],proj)

Вы также можете установить central_longitude=9.0 в своем TransverseMercator, это похоже на то, что указано в вашем шейп-файле.

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

Обновить

Кажется, что ваши границы также были установлены на основе только первого из municipalities:

In [34]: municipalities[0].bounds
Out[34]: (32458044.649189778, 5556418.748046352, 32465287.307457082, 5564153.5456742775)

Но другие элементы имеют другие границы. Вы можете сбросить лимиты до фактического чертежа на основе минимальных / максимальных значений границ всех municipalities.

bnd = np.array([i.bounds for i in municipalities])
x0,x1 = np.min(bnd[:,0]),np.max(bnd[:,2])
y0,y1 = np.min(bnd[:,1]),np.max(bnd[:,3])
ax.set_extent([x0,x1,y0,y1],proj)
person tsj    schedule 09.10.2016
comment
Спасибо за ваши усилия, @tsj. Я не совсем понимаю, почему границы шейп-файла, который я получил с помощью fiona, не дают мне точно наименьшего возможного квадрата, в котором может поместиться страна. Проекция (все еще) неверно указана в моем коде сверху? Есть ли ошибка в фионе или картопии? Хотя бы несогласованность? - person Jhonny; 10.10.2016
comment
Кроме того, просто подумайте, +/- 2e7 должно быть достаточно со стороны картографии, так как окружность Земли составляет около 40e6. Может быть, есть способ изменить ваш ввод так, чтобы он ссылался на -7.5e6 вместо 32.5e6, тогда вам не пришлось бы редактировать источник картопии. - person tsj; 10.10.2016
comment
Итак, метод геометрии фактически возвращает фигурный объект, поэтому я могу вызвать bounds, верно? И поскольку это возвращает (x0, y0, x1, y1), в то время как cartopy ожидает (x0,x1,y0,y1), мне нужно реорганизовать, как вы показываете во фрагменте кода. Понял, спасибо - работает как шарм. - person Jhonny; 13.10.2016