Matplotlib 底图海岸坐标

Question

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

Answer 1

从 drawcoastlines() 获取坐标的最佳方法是使用其 class 属性 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 公里。

Answer 2

这是另一种不依赖底图投影并提供原始 lon/lat 坐标的可能性。 advantage/disadvantage 是，大陆线未在地图边界处分割。

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()

Matplotlib 底图海岸坐标

Matplotlib Basemap Coastal Coordinates

python

matplotlib

matplotlib-basemap