Как узнать, находится ли координата внутри многоугольника в шейп-файле, который не работает?

У меня есть координаты (широта и долгота) в документе Excel. У меня есть шейп-файл, который содержит все фигуры различных провинций Канады. Я хотел бы иметь возможность создать новое поле в Excel, чтобы классифицировать разные координаты по разным провинциям Канады. Я попробовал приведенный ниже код, но он не работает.

import fiona
import shapely.geometry 

with fiona.open(r"D:\Users\Jonathan\Desktop\CRA-Project v2\Census Division\lcd_000b16a_e.shp") as fiona_collection:

shapefile_record = fiona_collection.next()

# Use Shapely to create the polygon
shape = shapely.geometry.asShape(shapefile_record['geometry'])
#print(shape)

point = shapely.geometry.Point(46.362914,-63.503809) # longitude, latitude

# Alternative: if point.within(shape)
if shape.contains(point):
    print("Found shape for point.")

Обновление 1: точка = shapely.geometry.Point (46.362914, -63.503809)

Многоугольник: Ссылка

Обновление решения. Я обновляю этот пост, потому что нашел решение и, надеюсь, оно кому-то поможет!

  1. Если вы используете шейп-файлы, с ними связана проекция (документация)
  2. Изучите документацию и определите, какую проекцию она использует. Например: в Канаде это конформная коника Ламберта.
  3. Преобразуйте ваши широты и долготы в эквивалент проекции
  4. Прокрутите шейп-файл, чтобы определить, находится ли новый эквивалент широты / долготы внутри многоугольника / мультиполигона. Вы можете сделать это с помощью Python!

python shapely
person Jonathan Lam    schedule 29.11.2018    source источник
comment
можешь поделиться тем, что не работает? (например: сообщение об ошибке)   -  person iamanigeeit    schedule 30.11.2018
comment
Указанные выше координаты находятся в пределах многоугольника, но я не получаю сообщение для печати. Таким образом, это не определение координат внутри многоугольника.   -  person Jonathan Lam    schedule 30.11.2018
comment
вы поменяли местами широту / долготу? см. gis.stackexchange.com/questions/84114/   -  person iamanigeeit    schedule 30.11.2018
comment
Цикл по всем полигонам очень затратен, если у нас много точек и много полигонов. Какая быстрая альтернатива?   -  person mah65    schedule 08.11.2020

Ответы (1)


arrow_upward
0
arrow_downward

С этим есть несколько проблем. Во-первых, ваш шейп-файл (мультиполигон) кажется в другой проекции, поэтому для работы в координатах долгота / широта необходимо ее преобразовать. Более того, даже после преобразования координата x является долготой, поэтому координаты точки, с которой вы проводите тестирование, должны быть поменяны местами. Пример может выглядеть так, как показано ниже. Предполагается, что входная проекция равна PCS_Lambert_Conformal_Conic (EPSG: 3347) - это указано в сопроводительном файле prj.

from functools import partial
import sys

import fiona
from shapely.geometry import Point, Polygon, asShape
from shapely.ops import transform
from shapely.wkt import loads
import pyproj

project = partial(
    pyproj.transform,
    pyproj.Proj(init='epsg:3347'),
    pyproj.Proj(init='epsg:4326'))

P = Point(-63.503809, 46.362914)
with fiona.open(sys.argv[1]) as F:
    for idx,feature in enumerate(F):
        G = transform(project, asShape(feature['geometry']))
        if G.contains(P):
            print(feature['properties'])
            break

Это производит:

OrderedDict([('CDUID', '1102'), ('CDNAME', 'Queens'), ('CDTYPE', 'CTY'), ('PRUID', '11'), ('PRNAME', 'Prince Edward Island / Île-du-Prince-Édouard')])

то есть он находит координаты в пределах Острова Принца Эдуарда.

person ewcz    schedule 30.11.2018