Matplotlib Basemap Coastal Coordinates

Question

Is there a way to query basemap to extract all coastal coordinates? Say user provides lat/lng and the function returns true/false if the coordinates are within 1km from the coast?

Answer 1

The best way to get the coordinates from drawcoastlines() is using its class attribute get_segments(). There is an example how you can get the distance from coast for a single point with longitude ans latitude in decimal degrees. You can adapt this function to use a unique map to calculate all points in a list. I hope it's help 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've got 0.7 Km.

Answer 2

This is another possibilities that doesn't rely on the basemap projection and gives raw lon/lat coordinates. An advantage/disadvantage is, that the continent lines are not split at the map boundaries.

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