2017-09-01 13 views
0

J'essaie de géoréférencer certaines données dans un .shp, en utilisant Geopandas.Trouver la distance à l'intérieur sur une arête de géométrie Geopandas, en fonction des points et du relèvement

Je voudrais trouver un point dans une arête de géométrie, compte tenu d'un relèvement et d'un point initial à l'intérieur de celle-ci.

Comment est-ce que je peux le faire d'une manière optimisée, puisque l'algorithme traitera beaucoup d'itérations (approximativement: dix par point)?

Image

Répondre

1

Cela dépend si vous voulez coller à LON/lat ou vous pouvez passer à système de coordonnées projetées comme UTM; et aussi ce que vous voulez dire exactement par bearing. Une voie possible en supposant que les coordonnées projetées et compass bearing (degré d'horloge du nord: 0-360) est de tracer une ligne dans la direction de votre roulement qui est assez long pour croiser le polygone et calculer la coordonnée de l'intersection.

disons que nous avons un GeoDataFrame contenant un Polygon d'un quartier à Berlin:

berlin 
# 
# bbox_east bbox_north bbox_south bbox_west geometry place_name 
# 0 13.429402 52.540407 52.504037 13.36586 POLYGON ((389163.2519209321 5821873.324153989,... Mitte, Berlin, Deutschland 

calculer x/y du centre de gravité de la géométrie (peut être un point quelconque, j'utilise ici barycentre pour représenter la idée):

x = berlin.centroid.geometry.item().x 
y = berlin.centroid.geometry.item().y 

fonction suivante peut calculer les coordonnées du nouveau point:

from shapely.geometry import LineString, LinearRing 

def find_point(polygon, x, y , bearing): 

    east, south, west, north = polygon.bounds 
    line_length = max(abs(east-west), abs(north-south)) * 2 

    new_x = x + (np.sin(np.deg2rad(bearing)) * line_length) 
    new_y = y + (np.cos(np.deg2rad(bearing)) * line_length) 

    l = LineString([[x,y], [new_x, new_y]]) 
    lr = LinearRing(polygon.exterior.coords) 
    intersections = lr.intersection(l) 

    return intersections.x, intersections.y 

essayer avec nos données d'entrée:

x_, y_ = find_point(berlin.geometry[0], x, y, 120) 

fig, ax = plt.subplots() 

berlin.plot(ax=ax) 
plt.scatter(x_, y_) 
plt.scatter(x,y) 

enter image description here

+0

Exactement ce que je cherchais. En effet, l'algorithme reçoit des données dans WGS84, puis les convertit en UTM. Le roulement est exactement ce que vous avez répondu! J'essaie d'apprendre Python à géoprocesser et faire des calculs pour les analyses de vent et d'onde, et parfois, il est difficile de trouver des informations en plus de la documentation de certaines bibliothèques. Merci beaucoup pour votre aide! –