Попытка создать обычную географическую сетку с помощью python — широта/долгота превышает ограничения

Я хочу получить все пары координат широты/долготы регулярной сетки в США. Я получил следующий код из более раннего вопроса:

import shapely.geometry
import pyproj

# https://stackoverflow.com/questions/40342355/how-can-i-generate-a-regular-geographic-grid-using-python
# Set up projections
p_ll = pyproj.Proj(init='epsg:4326') #4326
p_mt = pyproj.Proj(init='epsg:3857') # espg:3857 metric; same as EPSG:900913

# Create corners of rectangle to be transformed to a grid
nw = shapely.geometry.Point((24.0,-89))  #actually SW
print(nw)
se = shapely.geometry.Point((25.0,-65)) #actually NE

# 24.0,-89
# 50.0,-65

stepsize = 50000 # 5 km grid step size

# Project corners to target projection
s = pyproj.transform(p_ll, p_mt, nw.x, nw.y) # Transform NW point to 3857
e = pyproj.transform(p_ll, p_mt, se.x, se.y) # .. same for SE
print(s)
print(e)
# Iterate over 2D area
gridpoints = []
x = s[0]
while x < e[0]:
    y = s[1]
    while y < e[1]:
        p = shapely.geometry.Point(pyproj.transform(p_mt, p_ll, x, y))
        gridpoints.append(p)
        y += stepsize
    x += stepsize

# with open('testout.csv', 'wb') as of:
#     of.write('lon;lat\n')
for p in gridpoints:
    print(str(p.x)+"," +str(p.y))
    # of.write('{:f};{:f}\n'.format(p.x, p.y))

Источник: Как я могу создать обычную географическую сетку с помощью python ?

Код не работает, когда моя долгота превышает -89,9999 (от -90 до -180), что приводит к следующей ошибке:

MacBook-Pro:Desktop user$ python updated_calculate_points.py 
Traceback (most recent call last):
  File "updated_calculate_points.py", line 19, in <module>
    s = pyproj.transform(p_ll, p_mt, nw.x, nw.y) # Transform NW point to 3857
  File "/Users/user/anaconda3/lib/python3.6/site-packages/pyproj/__init__.py", line 504, in transform
    _proj._transform(p1,p2,inx,iny,inz,radians)
  File "_proj.pyx", line 361, in _proj._transform (_proj.c:3762)
RuntimeError: b'latitude or longitude exceeded limits'

Я не совсем понимаю, как работают эти прогнозы. Может ли кто-нибудь помочь?


person user3200950    schedule 09.12.2017    source источник
comment
Зачем вам нужна широта, превышающая +/- 90 градусов? (к северу от северного полюса или к югу от южного полюса) Это не допускается для большинства проекций.   -  person geocodezip    schedule 10.12.2017
comment
@geocodezip. Вы правы, ошибка связана с долготой между -90 и -180.   -  person user3200950    schedule 10.12.2017


Ответы (1)


Использование библиотеки h3 в сочетании с определенными географическими регионами (т. взято из Openstreetmap) — элегантное решение

person user3200950    schedule 23.08.2019