Matplotlib Graphics Cards Coastal Coordinates

Question

Is there a way to request a basemap to retrieve all coastal coordinates? Say the user provides lat/lng, and the function returns true/falseif the coordinates are within 1 km of the coast?

Answer 1

The best way to get the coordinates from drawcoastlines()is to use its class attribute get_segments(). There is an example of how you can get the distance from the coast to one point with longitude and latitude in decimal degrees. You can adapt this function to use a unique map to calculate all the points in the list. Hope this helps you.

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

Another way to get it:

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

I have 0.7 km.

Answer 2

This is another feature that does not depend on the projection of the base map and gives the original coordinates lon / lat. The advantage / disadvantage is that the lines of the continents are not divided at the borders of the map.

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