Matplotlib 底图沿海坐标答案

【问题标题】：Matplotlib Basemap Coastal CoordinatesMatplotlib 底图沿海坐标
【发布时间】：2015-12-27 12:18:42
【问题描述】：

有没有办法查询底图以提取所有海岸坐标？假设用户提供lat/lng，如果坐标在距离海岸1公里以内，函数返回true/false？

【问题讨论】：

有人有链接吗？为什么你会在没有任何解释的情况下对一个问题投反对票？

标签： python matplotlib matplotlib-basemap

【解决方案1】：

从drawcoastlines() 获取坐标的最佳方法是使用其类属性get_segments()。有一个例子，你可以得到一个点到海岸的距离，它的经度和纬度都是十进制度。您可以调整此函数以使用唯一地图来计算列表中的所有点。希望对你有帮助。

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np


def distance_from_coast(lon,lat,resolution='l',degree_in_km=111.12):
    plt.ioff()

    m = Basemap(projection='robin',lon_0=0,resolution=resolution)
    coast = m.drawcoastlines()

    coordinates = np.vstack(coast.get_segments())
    lons,lats = m(coordinates[:,0],coordinates[:,1],inverse=True)

    dists = np.sqrt((lons-lon)**2+(lats-lat)**2)

    if np.min(dists)*degree_in_km<1:
      return True
    else:
      return False

另一种获取方式：

from mpl_toolkits.basemap import Basemap
import matplotlib.pyplot as plt
import numpy as np
import os

def save_coastal_data(path,resolution='f'):

    m = Basemap(projection='robin',lon_0=0,resolution=resolution)

    coast = m.drawcoastlines()

    coordinates = np.vstack(coast.get_segments())
    lons,lats = m(coordinates[:,0],coordinates[:,1],inverse=True)

    D = {'lons':lons,'lats':lats}

    np.save(os.path.join(path,'coastal_basemap_data.npy'),D)

def distance_from_coast(lon,lat,fpath,degree_in_km=111.12):

    D = np.load(fpath).tolist()

    lons,lats = D['lons'],D['lats']

    dists = np.sqrt((lons-lon)**2+(lats-lat)**2)

    print np.min(dists)*degree_in_km

#Define path
path = 'path/to/directory'
#Run just one time to save the data. Will cost less time
save_coastal_data(path,resolution='h')  

distance_from_coast(-117.2547,32.8049,
os.path.join(path,'coastal_basemap_data.npy'))

我有 0.7 公里。

【讨论】：

这不是答案吗？
对不起，我刚从某个地方回来。看起来不错，我去看看
嗨，我有时间测试你的代码。这是我的坐标：-117.2547,32.8049，它返回 12.0341337387 公里。到目前为止这很好，但并不准确。在这里查看：gps-coordinates.net。知道如何使它更准确吗？我使用了 resolution='c'，它给了我 53.xx 公里？

【解决方案2】：

这是另一种不依赖底图投影并提供原始经度/纬度坐标的可能性。一个优点/缺点是，大陆线没有在地图边界处分割。

import matplotlib.pyplot as plt
from mpl_toolkits import basemap
import numpy as np
import os

def get_coastlines(npts_min=0):
    # open data and meta data files
    dirname_basemap = os.path.dirname(basemap.__file__)
    path_points = os.path.join(dirname_basemap, 'data', 'gshhs_c.dat')
    path_meta = os.path.join(dirname_basemap, 'data', 'gshhsmeta_c.dat')

    # read points for each segment that is specified in meta_file
    points_file = open(path_points, 'rb')
    meta_file = open(path_meta,'r')
    segments = []
    for line in meta_file:
        # kind=1 are continents, kind=2 are lakes
        kind, area, npts, lim_south, lim_north, startbyte, numbytes,\
        date_line_crossing = line.split()
        points_file.seek(int(startbyte))
        data = np.fromfile(points_file, '<f4', count = int(numbytes)/4)
        data = data.reshape(int(npts), 2)
        if npts_min < int(npts):
            segments.append(data)
    return segments


def main():
    segments = get_coastlines()
    fig, ax = plt.subplots(1, 1)
    for seg in segments:
        plt.plot(seg[:, 0], seg[:, 1])
    plt.show()


if __name__ == "__main__":
    main()

【讨论】：