geopandasmatplotlib-basemapshapelycartopypygmt

Plot milepost along coastline in python

I want to plot mileposts on a map every 100 miles along the coastline. An example is shown in the figure below: enter image description here

It is easy to plot the coastlines using Cartopy, but how to determine the locations of these mileposts and show them on a map? It would be better if the coastlines between the points were colored.

Solution

  • You can use Shapely's interpolate method to figure out the location of the mileposts. However, the tricky part of this is getting a singlepart linestring for the coastline. I messed around with a few coastline shapefiles that I downloaded, but due to the complexity and distance, getting a nice singlepart linestring was not a simple task. Therefore, I chose to digitize my own for this example (digitize.shp).

    import geopandas as gpd
import numpy as np
import matplotlib
import cartopy.crs as ccrs
import matplotlib.pyplot as plt

# hand digitized simplified version of coastline
cl_gdf = gpd.read_file("digitize.shp")


# project data from geographic crs so we have meters for interpolation below, vs degrees
cl_gdf = cl_gdf.to_crs(9311)

# get shapely ls from gdf
coastline = cl_gdf.iloc[0].geometry

interval = 160934 # approx 100 miles in meters
interval_arr = np.arange(0, coastline.length, interval )
points = [coastline.interpolate(d) for d in interval_arr] + [coastline.boundary[1]]

# create gdf from our list of interpolated Shapely points
points_gdf = gpd.GeoDataFrame(geometry=points, crs=9311)

# transform crs to wgs84 for plotting
points_gdf = points_gdf.to_crs(4326)

# add Lat and Long cols from geometry for plotting
points_gdf['Lat'] = (points_gdf['geometry'].apply(lambda geom: np.max([coord[1] for coord in geom.coords])))
points_gdf['Long'] = (points_gdf['geometry'].apply(lambda geom: np.max([coord[0] for coord in geom.coords])))

ax = plt.axes(projection=ccrs.PlateCarree())
ax.coastlines()

ax.set_extent([points_gdf['Long'].min() - 1,
               points_gdf['Long'].max() + 1,
               points_gdf['Lat'].min() - 1,
               points_gdf['Lat'].max() + 1], 
              crs=ccrs.PlateCarree())
              
# add our interpolated points to the Cartopy coastline
ax.scatter(points_gdf['Long'], points_gdf['Lat'], color='red', s=10)

fig = matplotlib.pyplot.gcf()
fig.set_size_inches(18.5, 10.5)

plt.show()

    enter image description here

    EDIT:

    Here is a replacement for "digitize.shp" for a full working example. You can use gpd's to_file() if you want to save cl_gdf to a shapefile.

    from shapely.geometry import LineString

ls_coords = [(-84.35698841221803, 29.95854398344339),
 (-84.05368623055452, 30.09279249008134),
 (-83.96252983715839, 30.0613020996354),
 (-83.7586709937452, 29.908822314318225),
 (-83.44376708928581, 29.68673219222582),
 (-83.40854757365547, 29.657727885236138),
 (-83.16946921461201, 29.28191493609843),
 (-82.7551219719023, 28.99684403311415),
 (-82.66893774541866, 28.698514018363166),
 (-82.68219685718537, 28.439961338912305),
 (-82.79489930720241, 28.158205213869703),
 (-82.63910474394356, 27.94937420354401),
 (-82.54629096157659, 27.641099854967987),
 (-82.63578996600191, 27.392491509342157),
 (-82.26122005859231, 26.72290636512327),
 (-81.52533935553987, 25.88095276793714),
 (-80.85243943337927, 25.168275510476434),
 (-80.7596256510123, 25.14175728694301),
 (-80.65355275687861, 25.150044231797207),
 (-80.59222936495757, 25.176562455330625),
 (-80.49278602670725, 25.206395456805726),
 (-80.65355275687864, 24.90143588617138),
 (-80.53587813994908, 25.03236961486765),
 (-80.23920551416893, 25.340643963443675),
 (-80.27069590461487, 25.353903075210386),
 (-80.37345402080688, 25.312468350939415),
 (-80.104957007531, 25.963822216479077),
 (-80.07180922811422, 26.91847826368225),
 (-80.66846925761622, 28.60238545805451),
 (-80.72150570468307, 28.86756769338872),
 (-81.19883372828465, 29.663114399391368),
 (-81.43749774008545, 30.392365546560455),
 (-81.47727507538558, 31.095098470196124),
 (-81.26512928711821, 31.406687596713834),
 (-81.09276083415097, 31.844238285015287),
 (-80.7612830399832, 32.175716079183054),
 (-80.48947124876561, 32.48067564981741),
 (-80.07180922811422, 32.61326676748452),
 (-79.32929896917841, 33.07070612343603),
 (-79.20996696327803, 33.22981546463656),
 (-78.9779325073606, 33.61432970587117),
 (-78.55364093082585, 33.8331050500219),
 (-77.93046267779044, 33.919289276505516),
 (-77.71831688952308, 34.28391485009006),
 (-77.45976421007221, 34.469542414824005),
 (-75.76259790393324, 35.20542311787645),
 (-75.76259790393324, 36.21974516802982),
 (-75.84215257453349, 36.58437074161438),
 (-76.17031559075959, 36.939051981373886),
 (-76.29959193048501, 36.96722759387815),
 (-76.2946197635725, 37.12136476816616),
 (-76.41726654741457, 37.16942904832049),
 (-76.50179338492735, 37.23738199612488),
 (-76.3377118768143, 37.5108511763133),
 (-76.72885567393226, 37.79923685723926),
 (-76.76200345334904, 37.878791527839525),
 (-76.31119365328087, 37.68653440722222),
 (-76.23163898268061, 37.89205063960623),
 (-76.35428576652268, 37.951716642556434),
 (-76.38411876799778, 37.95668880946896),
 (-76.50676555183985, 38.02298436830251),
 (-76.52665421948991, 38.05613214771928),
 (-76.6028941121485, 38.10253903890277),
 (-76.60952366803185, 38.12905726243619),
 (-76.85978940262852, 38.16717720876549),
 (-76.90785368278284, 38.19203804332807),
 (-76.94265885117046, 38.2102693220073),
 (-77.01724135485821, 38.313027438199306),
 (-77.01724135485821, 38.31799960511182),
 (-77.26916447842571, 38.3428604396744),
 (-77.30396964681333, 38.38263777497453),
 (-77.29568270195914, 38.52517322646668),
 (-77.23270192106726, 38.60804267500862),
 (-77.21612803135888, 38.63953306545456),
 (-77.17966547400042, 38.613014841921135),
 (-76.96254751882053, 38.4108133874788),
 (-76.81338251144503, 38.27987965878253),
 (-76.41063699153119, 38.311370049228465),
 (-76.53825594228582, 38.710800791200626),
 (-76.53494116434413, 39.2047027045106),
 (-76.0808165863343, 39.532865720736694),
 (-76.02446536132577, 39.370441601594486),
 (-76.06755747456758, 39.254424373635764),
 (-76.14048258928449, 39.101944588318595),
 (-76.1835747025263, 38.956094358884776),
 (-76.196833814293, 38.88316924416787),
 (-76.18688948046797, 38.76383723826747),
 (-76.17031559075959, 38.64782001030875),
 (-76.094075698101, 38.51854367058332),
 (-75.98468802602564, 38.33291610584937),
 (-75.89850379954201, 38.24341710142407),
 (-75.84546735247517, 38.190380654357234),
 (-75.77917179364162, 38.10419642787361),
 (-75.69961712304135, 38.011382645506636),
 (-75.65652500979955, 37.95171664255644),
 (-75.66978412156625, 37.91525408519799),
 (-75.92502202307544, 37.54068417778841),
 (-75.76591268187491, 37.481018174838205),
 (-75.67972845539128, 37.46112950718814),
 (-75.62337723038277, 37.48433295277989),
 (-75.58359989508264, 37.59703540279693),
 (-75.2189743214981, 38.051159980806766),
 (-75.05655020235588, 38.37600821909118),
 (-75.03666153470581, 38.44893333380809),
 (-75.39465755240701, 39.1450367015604),
 (-75.50404522448237, 39.39033026924455),
 (-75.54713733772417, 39.536180498678355),
 (-75.5438225597825, 39.60413344648275),
 (-75.50901739139488, 39.57595783397849),
 (-75.51896172521991, 39.48977360749487),
 (-75.50073044654069, 39.45331105013641),
 (-75.45100877741551, 39.41684849277796),
 (-75.3996297193195, 39.38701549130286),
 (-75.33333416048596, 39.34889554497357),
 (-75.2753255465066, 39.31409037658595),
 (-75.17588220825627, 39.27431304128582),
 (-75.07146670309342, 39.23453570598569),
 (-74.95710686410554, 39.19144359274388),
 (-74.8957834721845, 39.169897536122974),
 (-74.53778745448334, 39.270998263344154),
 (-74.29249388679919, 39.46325538396146),
 (-74.10023676618188, 39.89417651637956),
 (-74.02731165146497, 40.046656301696736),
 (-74.13338454559866, 40.32509764879766),
 (-74.23945743973235, 40.4974661017649),
 (-74.07371854264846, 40.550502548831744),
 (-73.80853630731424, 40.56376166059845),
 (-73.60302007493023, 40.550502548831744),
 (-73.4439107337297, 40.61016855178194),
 (-73.39750384254621, 40.70298233414891),
 (-73.41739251019628, 40.88198034299951),
 (-73.42070728813798, 40.935016790066356),
 (-73.62953829846366, 40.8985542327079),
 (-73.77207374995581, 40.83225867387435),
 (-73.7588146381891, 40.931702012124674),
 (-73.53672451609668, 41.031145350375006),
 (-73.15552505280375, 41.1504773562754),
 (-72.96658271012812, 41.17368080186715),
 (-71.70033753640726, 41.35930836660109),
 (-71.04401150395508, 41.50515859603491),
 (-70.70590415390394, 41.65100882546872),
 (-70.68601548625388, 41.670897493118794),
 (-70.67938593037053, 41.51178815191826),
 (-70.61309037153697, 41.53830637545168),
 (-70.5070174774033, 41.76371127548577),
 (-70.54016525682006, 41.86978416961945),
 (-70.6528677068371, 42.00237528728656),
 (-70.77882926862085, 42.247668854970726),
 (-70.93130905393804, 42.41340775205461),
 (-70.91142038628796, 42.5924057609052),
 (-70.86501349510448, 42.65870131973875),
 (-70.54679481270341, 43.32165690807428),
 (-70.19542835088558, 43.63987559047534),
 (-69.74461855081742, 43.81887359932593),
 (-69.53910231843341, 43.885169158159485),
 (-68.96233095658148, 44.289572067044155),
 (-68.81648072714766, 44.44868140824468),
 (-68.41207781826299, 44.48845874354482),
 (-67.91486112701133, 44.42216318471126),
 (-67.61653111226035, 44.52823607884495),
 (-67.09942575335863, 44.70723408769554),
 (-67.02650063864172, 44.760270534762384)]

ls = LineString(ls_coords)
cl_gdf = gpd.GeoDataFrame(geometry=[ls], crs=4326)