I want to plot mileposts on a map every 100 miles along the coastline. An example is shown in the figure below:
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.
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()
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)
