判断经纬度是否位于陆地
最近，师姐问我如何去除落在陆地上的经纬度点。巧的是，我之前正好处理过类似的问题。解决这个问题的关键在于判断给定的经纬度是否位于陆地上，从而筛选出符合要求的数据。
TC生成点图
首先我们读取数据，该数据是使用热带气旋(TC)历史最佳路径数据(IBTrACS)获得，我们记录了印度洋上TC生成位置以及类别，保存为csv格式(文件名为Indian_Tropical_Cyclogenesis.csv)，里面的内容如下：
为了方便，这里我们只展示部分数据内容。
绘制结果如下：
注意：
图中用红色圈标出了位于陆地上的TC生成点，但我们需要筛选出全部位于海洋上的TC生成点。
绘图代码
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import pandas as pd
import matplotlib.patches as mpatches
import os
f = pd.read_csv(r'Indian_Tropical_Cyclogenesis.csv')
def get_style(grade):
color_map = {
'TD': 'lime',
'TS': 'mediumblue',
'STS': 'darkorange',
'TY': 'crimson',
'STY': 'fuchsia'
}
marker_map = {
'TD': 's',
'TS': '^',
'STS': '>',
'TY': 'D',
'STY': '+'
}
color = color_map.get(grade, 'purple')
marker = marker_map.get(grade, 'v')
return color, marker
fig = plt.figure(figsize=(10,8),dpi=300)
proj = ccrs.PlateCarree()
region = [45, 105, 0, 32]
ax = plt.axes(projection=proj)
ax.set_extent(region,crs = proj)
land = cfeature.NaturalEarthFeature('physical', 'land', '50m')
ax.add_feature(land, color='lightgray', zorder=-1)
for i in range(len(f['year'])):
tc_lon = f['lon'][i]
tc_lat = f['lat'][i]
grade = f['grade'][i]
color, marker = get_style(grade)
fi_color = str(color)
fi_marker = str(marker)
point = sgeom.Point(tc_lon, tc_lat)
sc = ax.scatter(tc_lon,tc_lat,c=fi_color,marker=fi_marker,s=30)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.LAND, color='lightgray', zorder = -1)
ax.add_feature(cfeature.OCEAN, color='whitesmoke')
ax.set_xticks([55,65,75,85,95])
ax.set_yticks([0,10,20, 30])
ax.set_xticklabels(['55°E','65°E','75°E','85°E','95°E'],fontsize=15)
ax.set_yticklabels(['0°','10°N','20°N','30°N'],fontsize=15)
ax.set_title('1991-2020 TC source', fontsize=15, loc='left', fontweight = "bold")
colors = ['lime', 'mediumblue', 'darkorange','crimson', 'fuchsia', 'purple']
texts = ['TD', 'TS', 'STS', 'TY', 'STY', 'SSTY']
patches = [mpatches.Patch(color=colors[i],label=texts[i]) for i in range(len(colors))]
for i in range(len(texts)):
color, markers = get_style(texts[i])
ax.scatter([], [], c = colors[i], marker = str(markers), label = texts[i])
plt.legend(fontsize =10)
plt.show()
解决方案
通过将经纬度点转换为几何点对象，并检查这些点是否位于预先创建的陆地多边形内，从而判断它们是否位于陆地上。
Core
land_polygons = []
for poly in land.geometries():
if isinstance(poly, sgeom.Polygon):
land_polygons.append(sgeom.Polygon(list(poly.exterior.coords)))
elif isinstance(poly, sgeom.MultiPolygon):
for subpoly in poly.geoms:
land_polygons.append(sgeom.Polygon(list(subpoly.exterior.coords)))
for i in range(len(f['year'])):
tc_lon = f['lon'][i]
tc_lat = f['lat'][i]
grade = f['grade'][i]
color, marker = get_style(grade)
fi_color = str(color)
fi_marker = str(marker)
point = sgeom.Point(tc_lon, tc_lat)
if any([poly.contains(point) for poly in land_polygons]):
continue
sc = ax.scatter(tc_lon,tc_lat,c=fi_color,marker=fi_marker,s=30)
完整代码
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import cartopy.feature as cfeature
import pandas as pd
import shapely.geometry as sgeom
import matplotlib.patches as mpatches
import os
f = pd.read_csv(r'Indian_Tropical_Cyclogenesis.csv')
def get_style(grade):
color_map = {
'TD': 'lime',
'TS': 'mediumblue',
'STS': 'darkorange',
'TY': 'crimson',
'STY': 'fuchsia'
}
marker_map = {
'TD': 's',
'TS': '^',
'STS': '>',
'TY': 'D',
'STY': '+'
}
color = color_map.get(grade, 'purple')
marker = marker_map.get(grade, 'v')
return color, marker
fig = plt.figure(figsize=(10,8),dpi=300)
proj = ccrs.PlateCarree()
region = [45, 105, 0, 32]
ax = plt.axes(projection=proj)
ax.set_extent(region,crs = proj)
land = cfeature.NaturalEarthFeature('physical', 'land', '50m')
ax.add_feature(land, color='lightgray', zorder=-1)
land_polygons = []
for poly in land.geometries():
if isinstance(poly, sgeom.Polygon):
land_polygons.append(sgeom.Polygon(list(poly.exterior.coords)))
elif isinstance(poly, sgeom.MultiPolygon):
for subpoly in poly.geoms:
land_polygons.append(sgeom.Polygon(list(subpoly.exterior.coords)))
for i in range(len(f['year'])):
tc_lon = f['lon'][i]
tc_lat = f['lat'][i]
grade = f['grade'][i]
color, marker = get_style(grade)
fi_color = str(color)
fi_marker = str(marker)
point = sgeom.Point(tc_lon, tc_lat)
if any([poly.contains(point) for poly in land_polygons]):
continue
sc = ax.scatter(tc_lon,tc_lat,c=fi_color,marker=fi_marker,s=30)
ax.add_feature(cfeature.COASTLINE.with_scale('50m'))
ax.add_feature(cfeature.LAND, color='lightgray', zorder = -1)
ax.add_feature(cfeature.OCEAN, color='whitesmoke')
ax.set_xticks([55,65,75,85,95])
ax.set_yticks([0,10,20, 30])
ax.set_xticklabels(['55°E','65°E','75°E','85°E','95°E'],fontsize=15)
ax.set_yticklabels(['0°','10°N','20°N','30°N'],fontsize=15)
ax.set_title('1991-2020 TC source', fontsize=15, loc='left', fontweight = "bold")
colors = ['lime', 'mediumblue', 'darkorange','crimson', 'fuchsia', 'purple']
texts = ['TD', 'TS', 'STS', 'TY', 'STY', 'SSTY']
patches = [mpatches.Patch(color=colors[i],label=texts[i]) for i in range(len(colors))]
for i in range(len(texts)):
color, markers = get_style(texts[i])
ax.scatter([], [], c = colors[i], marker = str(markers), label = texts[i])
plt.legend(fontsize =10)
plt.show()
Result:
往期回顾
Python | MJO | 位相图 Python | 半球间经度转换 Python | 大气科学 | 偏相关 Python | 气象绘图 | 台风降水 Python | 解决 cmaps 库报错的问题 Python | 批量下载 NCEP2 再分析数据 Python | 北大西洋涛动 | NAO 指数 | EOF
