Как узнать, не работает ли координата внутри многоугольника в файле формы? - PullRequest
Новая
       6

Как узнать, не работает ли координата внутри многоугольника в файле формы?

0 голосов
/ 30 ноября 2018

У меня есть координаты (широта и долгота) в документе 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!

1 Ответ

0 голосов
/ 30 ноября 2018

Есть несколько проблем с этим.Во-первых, ваш шейп-файл (мультиполигон), похоже, находится в другой проекции, поэтому для работы в координатах долгота / широта необходимо преобразовать его.Более того, даже после преобразования координата 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')])

, т. Е.он находит координаты на острове Принца Эдуарда.

...