检查陆地上或海上的LAT／Lon点
global-land-mask的Python项目详细描述
全球陆地面罩
检查地球上任何一点的纬度/经度点是否在陆地上。
说明
这个python模块，global land mask，包含用于检查lat/lon点是在陆地上还是在海上的脚本。为了做到这一点，我们使用地球数据集，它以1km的分辨率对整个地球进行采样。然后，我们简单地从这个立面图中提取所有“无效”值并保存到文件中。
全局掩码是形状的（2160043200），在不压缩的情况下保存时约等于980mb。使用numpy savez_compressed可以将这些数据压缩到2.5 MB，从而使包非常紧凑。
全球数据集的原始高程数据可以从 https://www.ngdc.noaa.gov/mgg/topo/gltiles.html 无需下载此数据即可使用全局土地掩码。但是，通过下载，可以使用提供的函数构建全局高程数据集。
此包提供globe.is_land（），它是basemap.is_land（）的一个属性。对于10000个数据点，globe.is_land的速度大约是basemap.is_land的6000倍。
简单示例
下面是一个简单的例子，显示了使用global_land_mask检查纬度/经度点是否在陆地上。
fromglobal_land_maskimportglobeimportnumpyasnp# Check if a point is on land:lat=40lon=-120is_on_land=globe.is_land(lat,lon)print('lat={}, lon={} is on land: {}'.format(lat,lon,is_on_land))# lat=40, lon=-120 is on land: True# Check if several points are in the oceanlat=40lon=np.linspace(-150,-110,3)is_in_ocean=globe.is_ocean(lat,lon)print('lat={}, lon={} is in ocean: {}'.format(lat,lon,is_in_ocean))# lat=40, lon=[-150. -130. -110.] is in ocean: [ True True False]
速度测试
比较global_land_mask和basemap的性能。
fromglobal_land_maskimportglobefrommpl_toolkits.basemapimportBasemapimportmatplotlibmatplotlib.use('TkAgg')importnumpyasnpimporttime# Lat/lon points to getlat=np.linspace(-20,50,100)lon=np.linspace(-130,-70,100)# Make a gridlon_grid,lat_grid=np.meshgrid(lon,lat)# Get whether the points are on land using globe.is_landstart_time=time.time()globe_land_mask=globe.is_land(lat_grid,lon_grid)globe_run_time=time.time()-start_timeprint('Time to run globe.is_land(): {}'.format(globe_run_time))# Get whether the points are on land using Basemap.is_land# bm = Basemap(projection='cyl',resolution='i')bm=Basemap(projection='cyl',llcrnrlat=-60,urcrnrlat=90, \ llcrnrlon=-180,urcrnrlon=180,resolution='c')f=np.vectorize(bm.is_land)start_time=time.time()xpt,ypt=bm(lon_grid,lat_grid)basemap_land_mask=f(xpt,ypt)basemap_run_time=time.time()-start_timeprint('Time to run Basemap.is_land(): {}'.format(basemap_run_time))print('Speed up: {}'.format(basemap_run_time/globe_run_time))# Check agreement (note there is a different treatment for lakesfraction_agreed=np.sum(globe_land_mask==basemap_land_mask)/len(globe_land_mask.flatten())print('Fraction agreeing: {}'.format(fraction_agreed))
我们上方的地图示例
尝试运行
fromglobal_land_maskimportglobefrommpl_toolkits.basemapimportBasemapimportnumpyasnpimportmatplotlibmatplotlib.use('TkAgg')importmatplotlib.pyplotasplt# Lat/lon points to getlat=np.linspace(-20,90,1000)lon=np.linspace(-130,-60,1002)# Make a gridlon_grid,lat_grid=np.meshgrid(lon,lat)# Get whether the points are on land.z=globe.is_land(lat_grid,lon_grid)# Set up mapfig=plt.figure(0,figsize=(5.5,4.5))plt.clf()ax=fig.add_axes([0.1,0.1,0.8,0.8])m=Basemap(llcrnrlon=-119,llcrnrlat=22,urcrnrlon=-64,urcrnrlat=49,projection='lcc',lat_1=33,lat_2=45,lon_0=-95,area_thresh=200,resolution='i')m.drawstates(linewidth=0.2)m.drawcoastlines(linewidth=0.2)m.drawcountries(linewidth=0.2)cs=m.contourf(lon_grid,lat_grid,z,levels=[-0.5,0.5,1.5],cmap="jet",latlon=True)plt.show()plt.savefig('example_plot_globe_map_us.png',bbox_inches='tight',dpi=400)
要为美国创建二进制掩码，如下图所示：
请注意，湖泊被列为“陆地上的”，分辨率并不完美，但它足够好的许多目的！
欢迎加入QQ群-->： 979659372
