Пересечение двух ортодромов с использованием картопии - PullRequest
Новая
       23

Пересечение двух ортодромов с использованием картопии

0 голосов
/ 17 января 2020

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

lat1a = 49.9
lon1a = 22.5
lat1b = -20.2
lon1b = 11.99

lat2a = -51.5
lon2a = -68
lat2b = -8.04
lon2b = 14.6

land_50m = cfeature.NaturalEarthFeature('physical', 'land', '50m',
                                        edgecolor='face',
                                        facecolor=cfeature.COLORS['land'])

ax = plt.axes(projection=ccrs.Robinson())
ax.set_extent([-90, 50, -60, 70])
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.BORDERS, linewidth=0.2)


ax.plot(22.5, 49.2, marker='s', color='black', markersize=7, transform=ccrs.Geodetic())
ax.plot(-69.3, -51.5, marker='s', color='black', markersize=7, transform=ccrs.Geodetic())

ax.plot([lon2b, lon2a], [lat2b, lat2a], color='b', linestyle='--', transform=ccrs.Geodetic())
ax.plot(lon2b, lat2b, marker='o', color='red', markersize=5, transform=ccrs.Geodetic())
ax.plot([lon1b, lon1a], [lat1b, lat1a], color='b', linestyle='--', transform=ccrs.Geodetic())
ax.plot(lon1b, lat1b, marker='o', color='black', markersize=5, transform=ccrs.Geodetic())
plt.show()

и результат

enter image description here

I ' Мы думали о получении списка координат, которые находятся внутри синих пунктирных линий, и просто ищем повторение в обоих списках, чтобы я мог различить guish пересечение, но я застрял в этой части. Может быть, у вас есть другие идеи?

1 Ответ

1 голос
/ 18 января 2020

Чтобы получить то, что вы хотите (точка пересечения по Cartopy + shapely), требуется некоторое сложное кодирование, как вы можете видеть в предоставленном коде. Если вам просто нужно list of coordinates ..., просто посмотрите на bpath[0].vertices.

Рабочий код: 

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
from shapely.geometry import LineString

class LowerThresholdRobinson(ccrs.Robinson):
    # credits: https://stackoverflow.com/questions/40270990/cartopy-higher-resolution-for-great-circle-distance-line
    # use this to get finer geodesic line on Robinson projection
    @property
    def threshold(self):
        return 1e3

def get_Xpoint(b_handle, g_handle):
    """ Find (x,y) of the point of intersection """
    bpath = b_handle[0]._get_transformed_path().get_transformed_path_and_affine()
    gpath = g_handle[0]._get_transformed_path().get_transformed_path_and_affine()
    # Use vertices of the geodesic lines to  ...
    # create shapely LineString geometries
    bLs = LineString( bpath[0].vertices )
    gLs = LineString( gpath[0].vertices )

    answer = (None, None)

    # Find point of intersection between 2 LineString 
    if bLs.intersects(gLs): # returns True/False
        ans = bLs.intersection(gLs)
        answer = (ans.x, ans.y)

    return answer

# Begin:- Main code
lat1a = 49.9
lon1a = 22.5
lat1b = -20.2
lon1b = 11.99

lat2a = -51.5
lon2a = -68
lat2b = -8.04
lon2b = 14.6

# prep projections to get better solution
myProj = LowerThresholdRobinson()  # use this instead of `ccrs.Robinson()`
geodProj = ccrs.Geodetic()

fig = plt.figure(figsize=[6, 6])
ax = fig.add_subplot(1, 1, 1, projection=myProj)

#ax.set_extent([-90, 50, -60, 70])
ax.set_extent([0, 25, -22, 5])  #zoom-in to the intersection point
ax.add_feature(cfeature.COASTLINE, edgecolor='gray', facecolor='lightyellow')
ax.add_feature(cfeature.BORDERS, edgecolor='black', linewidth=0.2)

ax.plot(22.5, 49.2, marker='s', color='black', markersize=7, transform=geodProj)
ax.plot(-69.3, -51.5, marker='s', color='black', markersize=7, transform=geodProj)

b_handle = ax.plot([lon2b, lon2a], [lat2b, lat2a], color='b', linestyle='--', transform=geodProj)
ax.plot(lon2b, lat2b, marker='o', color='red', markersize=5, transform=geodProj)

g_handle = ax.plot([lon1b, lon1a], [lat1b, lat1a], color='g', linestyle='--', transform=geodProj)
ax.plot(lon1b, lat1b, marker='o', color='black', markersize=5, transform=geodProj)

# This gets and plots the point of intersection (as big brown + symbol)
xy = get_Xpoint(b_handle, g_handle)
ax.plot( xy[0], xy[1], marker='+', color='brown', markersize=60)

# This converts/prints xy to (long,lat)
print(geodProj.transform_point(xy[0], xy[1], myProj))

plt.show()

Выход:

enter image description here

...