This is a deep dive with a lot of repeative plotting to demonstrate how features are added to a graphical object and the impact each one makes to the entire product. The goal with this stragegy is that it will let us slowly explore how the more advanced graphics in Python work, albiet at the cost of a much longer Python notebook.
In this "Let's Play" session, we will explore the aspects of plotting graphics beyond the simple way that we've been doing it through the semester.
One reason for this is that when one goes to the "online 'help,'" the examples that are often provided for "how to make a plot" for example tends to expect a user to jump in and intuit what is happening.
Add the more abstract concepts of looking at map projections and how they are managed in any kind of software and speaking for me, personally, this is an engraved invitation to rage quit, which I did -- more than once -- when transitioning to Python. Therefore, in this session we are going to slow-roll the process of working with plots that have maps on them.
This tutorial is NOT about plotting meteorological or other geosiatial gridded data. Here we are ONLY going to do some rudementary plotting of maps with nothing more than a few points and lines on top of them.
For this exercize to get started it's a good idea to have these skills under your belt.
To begin we'll load our libraries with two previous old standards:
########################################################
#
# Python Libraries
#
#
# Numpy for arrays and basic math
#
import numpy as np
#
# MatPlotLib's Liraries
#
# Pyplot for basic plotting
import matplotlib.pyplot as plt
#
########################################################
Some of you may need to work with geospatial data. We actually can use the above example to show us the way.
We had declared a new "projection" for our polar plot called, "polar."
Now, imagine if you could now plot data in a different cooridate framework, such as latitude and longitude.
This will require us to load a new package, Cartopy. Which contains information as well as the means to plot basic map features like continents and countries.
Depending on what your application is you may need to install more than one package.
I recommend that it be done this way if you are using anaconda or miniconda.
!conda install shapely cartopyfrom inside Jupyter. But you will also need to restart your Jupyter Kernel for it to take hold
You can then load the suplemental libraries. You should only need the cartopy library to get this going, but you will need to install it much like when you load libaries from scipy.
Cartopy.crs: Cartopy's list of map projections. It's job is to connect the x and y position (we call them eastings and northings) to latitude and longitude. CRS stands for "Coordinate Reference System."
Cartopy.features: The features interface in Cartopy draws a set of pre-prepared geospatial features onto the Axes to where the projection (taken from Cartopy.crs) has been defined.
Implicitly loaded with Cartopy will also be the the Shapely Library.
Shapely's job is to manipilate and transform objects in a 2-d cartesian plan. Cartesian planes in this context are just "maps." For example, a map will have polygons, lines, etc corresponding to locatoins, borders, coastlines, etc, and later on as we explore, 2-D and 3-D grids. Typically these will be in a given projection such as Polar Stereographic, or Lambert Conformal Conic, or in terms of latitude and longitude. Shapely bends and warps these fields to fit a projection or to convert projections into latitude and longitude.
Note that Shapely is about to go through major new release. At the time of me writing this, the current version is 1.8. The next expected version, 2.0, will have a few major changes from the 1.* editions of Shapely. Therefore, when running things that will leverages these changes, even if you are not directly using the features, your notebook session will be dripping with so much pink you'll think Strawberry Shortcake killed the Klingon Chancelor with her teeth and claws alone. I'm trying to limit my use here of Cartopy and Shapely to the bare minium of what we need to plot in this exercise.
########################################################
#
# Cartopy Python Library
#
import cartopy.crs as ccrs
import cartopy.feature as cfeature
#
########################################################
To start let's drop an empty Axis with a projection onto our workspace.
For this example, I'm starting with a global projection called the "Interrupted Goode Homolosine" projection, a projection you may often see in atlases that resembles a peeled origins.
########################################################
#
# Cartopy Example: Laying down our plotting Canvas and Axes
#
# (1) Open a new default (single panel) figure canvas
# (2) Reset the axis at (0,0) to be in Robinson Projection Coordinates
# (3) Set a title for the plot
#
fig, ax = plt.subplots()
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.InterruptedGoodeHomolosine())
ax.set_title( "Interrupted Goode Homolosine")
#
########################################################
Text(0.5, 1.0, 'Interrupted Goode Homolosine')
Here apart from showing the outlines of the classic Goodd Homolosine "orange peel." That's ok. Because we really need to decide what the map is going to be of.
Will it just have continents? Will we want to zoom in? What features do we want to highlight borders, terrain height, something else?
But for now we are going to do a very simple line map. To do so, we will dig into the "features" library in Cartopy and just grab the coastlines.
To add a geographical feature use this example which leverages a command designed to use with Matplotlib, add.feature(). The link with this command includes information on some basic map global scale maps as well as US state boundaries.
(I am also going to change projections to the Robinson Projection which is a simple global coverage without having to "rip" the Earth's surface like we do with the Goode projection. (I'll be showing other projections as we play.)
########################################################
#
# Cartopy Example: Laying down our plotting Canvas and Axes
#
# (1) Open a new default (single pannel) figure canvas
# (2) Reset the axis at (0,0) to be in Robinson Projection Coordinates
# (3) Set a title for the plot
# (4) Render continental boundaries.
#
fig, ax = plt.subplots()
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.Robinson())
ax.add_feature(feature = cfeature.COASTLINE)
ax.set_title("Robinson Projection")
#
########################################################
Text(0.5, 1.0, 'Robinson Projection')
/opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:825: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:877: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:944: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0:
A map is nice... but at some point we want to put our own data on it.
Let's play with known places on the globe
To do this here is where it gets a little tricky. For example, plotting data in latitude and longitude would be intuitive. But in reality you are plotting in an x-y map space in which rather than latitude and longitude, you are plotting in "eastings" and "northings."
The good news is that your base background can use one projection scheme while you can plot data in another map projetion. Notice the syntax in the examples below. We are adding a "new" argument to the plot command: "transform," which you may remember can also be used in plotting in simple polar coordinates with r, and $\theta$.
########################################################
#
# Cartopy Example: Laying down our plotting Canvas and Axes
#
# (0) Get latitude and longitudes for various locations
#
# (1) Open a new default (single panel) figure canvas
# (2) Reset the axis at (0,0) to be in Robinson Projection Coordinates
# (3) Set a title for the plot
# (4) Render continental boundaries.
# (5) Plot The points in sequence.
# (6) Add Labels.
#
# Add Coordinates of Plotted Points (units are in decimal degrees north and east)
#
#
# There is a MUCH better way to do this particular plot in Pandas
# using a Data Frame and I'll add this at the very bottom.
#
KRAP = [-103.0605, 44.0384] # Rapid City, SD
KMSP = [ -93.2223, 44.8848] # Minneapolis, MN
KSLC = [-111.9791, 40.7899] # Salt Lake City, UT
KDEN = [-104.6737, 39.8561] # Denver, CO
KORD = [ -87.9090, 41.9803] # Chicago (O'Hare), IL
KDFW = [ -97.0403, 32.8998] # Dallas-Ft Worth, TX
YSSY = [ 151.2093, -33.8688] # Sydney, AU
SCEL = [ -70.7944, -33.3898] # Santiago, CL
locations = ["Rapid City",
"Minneapolis",
"Salt Lake City",
"Chicago",
"Denver",
"Dallas-Fort Worth",
"Sydney",
"Santiago"]
#
# Add the plot
#
fig, ax = plt.subplots(figsize=[11,8.5])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.InterruptedGoodeHomolosine())
ax.add_feature(feature = cfeature.COASTLINE)
ax.plot(KRAP[0], KRAP[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KMSP[0], KMSP[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KSLC[0], KSLC[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KDEN[0], KDEN[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KORD[0], KORD[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KDFW[0], KDFW[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(YSSY[0], YSSY[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(SCEL[0], SCEL[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.legend(locations)
ax.set_title("Interrupted Goode Homolosine")
plt.show()
#
########################################################
If you have heard about the "flat earther" craze, they project their world in what is often a polar projection.
Fun fact... the Soviets had this as their go-to projection as well which resulted in a very different worldview (literally) than United States who often uses a Mercator Projection, Robinson- or Goode-style projection when plotting things globaly. This comes as a surprise to a lot of people but in Meteorology, the Mercator projection is often a prefered projection in Numerical Weather Prediction (NWP) when working in the Tropics since a Mercator Projection is a "conformal" projection, meaning that it preserves angles -- i.e., vectors.
This leads to some genuine bemusent when they ahve to explain ceratin aircraft routes.
Let's repeat the previous graph in the Polar Stereographic Projection, which, BTW, is popilar as you move poleward in NWP.
Also from above, you can see that if you chose to point a set area of points, the total map extent will "clip" to the fit the data available. We can fix that with the axes.set_global() modifier (which acutally sits under Cartopy's library.
We are also going to add grid lines.
########################################################
#
# Cartopy Example: Laying down our plotting Canvas and Axes
#. (Cold War / Flat Earth Edition)
#
# (0) Get latitude and longitudes for various locations
#
# (We'll dispense with it here and use the ones from above)
# (1) Open a new default (single panel) figure canvas
# (2) Reset the axis at (0,0) to be in a Polar Stereographic Projection
# (3) Set a title for the plot
# (4) Render continental boundaries.
# (5) Plot The points in sequence.
# (6) Add Labels
# (7) Add Grid Lines
# (8) Set for a full global view.
fig, ax = plt.subplots(figsize=[11,8.5])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.NorthPolarStereo())
ax.add_feature(feature = cfeature.COASTLINE)
ax.plot(KRAP[0], KRAP[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KMSP[0], KMSP[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KSLC[0], KSLC[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KDEN[0], KDEN[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KORD[0], KORD[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(KDFW[0], KDFW[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(YSSY[0], YSSY[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.plot(SCEL[0], SCEL[1],
linestyle = "None",
marker = 'x',
transform = ccrs.PlateCarree())
ax.legend(locations)
ax.set_title("Polar Stereographic Projection")
ax.gridlines()
ax.set_global()
plt.show()
#
########################################################
/opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings.extend(multi_line_string) /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings.extend(multi_line_string)
Looking at the above problem, you can start to see the Sydney-to-Santiago flight from the Flat Earther's perspective. Traveling from Sydney to Santiago SHOULD be imposible given the lenght of travel needed in their worldview. Indeed, they believe that Qantus Flights 27 & 28, LATAM Flights 802 & 803 are myths and similar long hauls in the Southern Hemisphere are hoaxes. Rather, they argue that the real flights operate, for example as SYD->LAX->SCL -- which in reality reflects the economics of air travel given that the market for direct flights between Chile & Australia aren't as profitable (and thus, as common) as ~car~plane-pooling between more frequent trips like the daily trips passing through more global hubs like LAX.
So let's map out a few routes.
When we do this, you will want to have our line-of-sight paths using a different projection: Specifcally the Geodedic Projection. This will honor a spherical earth regardless of the projection you are using and provide a more credible "as the crow flies" distance.
########################################################
#
# Cartopy Example: Plotting Impossible Flight Routes!
#
# (0) Create our Flight Paths
#
# (1) Open a new default (single pannel) figure canvas
# (2) Reset the axis at (0,0) to be in Polar Stereographic Projection Coordinates
# (3) Set a title for the plot
# (4) Render continental boundaries.
# (5) Plot our flight paths using a Geodedic rendering!
# (6) Add Labels.
# (7) Add Grid Lines
#
# Bring out our flight paths!@
#
DL_2599 = np.array([KRAP,KMSP])
DL_3710 = np.array([KRAP,KSLC])
UA_4743 = np.array([KRAP,KDEN])
UA_5427 = np.array([KRAP,KORD])
AA_5765 = np.array([KRAP,KDFW])
QA___27 = np.array([YSSY,SCEL])
flights = ["DL_2599",
"DL_3710",
"UA_4743",
"UA_5427",
"AA_5765",
"QA___27"]
fig, ax = plt.subplots(figsize=[11,8.5])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.NorthPolarStereo())
ax.add_feature(feature = cfeature.COASTLINE)
ax.plot(DL_2599[:,0], DL_2599[:,1], transform = ccrs.Geodetic())
ax.plot(DL_3710[:,0], DL_3710[:,1], transform = ccrs.Geodetic())
ax.plot(UA_4743[:,0], UA_4743[:,1], transform = ccrs.Geodetic())
ax.plot(UA_5427[:,0], UA_5427[:,1], transform = ccrs.Geodetic())
ax.plot(AA_5765[:,0], AA_5765[:,1], transform = ccrs.Geodetic())
ax.plot(QA___27[:,0], QA___27[:,1], transform = ccrs.Geodetic())
ax.legend(flights)
ax.set_title("Polar Stereographic Projection")
ax.gridlines()
plt.show()
#
########################################################
/opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:825: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:877: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:944: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings.extend(multi_line_string) /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings.extend(multi_line_string)
Checkmate Flat Earthers!
While that last flight looks like a helluva miles run, remember that they're in the southern hemisphere. So let's redo those flights keeping our domestic flights domestic and that one big one using a more approrate projection such as the South Polar Stereographic Projection!
Also since the path will still be long and the projected straight line (on a sphere) looks like a curve in this projection, a careful look at the QA 27 flight path looks like a collection of line segments rather than an smooth curve. This is because the rendering of lines in projected space are based on a default map projection's distance threshold length -- the longest line it will draw. We will want small line segment lengths as we draw our circular paths.
To do this, the second plot below illustrates how that can be done: You must "clone" the projection that your MAP is on into a new variable. Then you can manipulate various parameters of the projection -- for us that is just the "threshold" attribute. The example is here and the second plot shows it in action.
hirezSouthPS = ccrs.SouthPolarStereo()
hirezSouthPS._threshold = hirezSouthPS._threshold/100.########################################################
#
# Cartopy Example: Plotting Domestic Flight Routes!
#
# (0) Create our Flight Paths
#
# (1) Open a new default (single pannel) figure canvas
# (2) Reset the axis at (0,0) to be in Polar Stereographic Projection Coordinates
# (3) Set a title for the plot
# (4) Render continental and US State boundaries.
# (5) Plot our flight paths using a Geodedic rendering!
# (6) Add Labels.
# (7) Add Grid Lines
#
# Bring out our flight paths!@
#
DL_2599 = np.array([KRAP,KMSP])
DL_3710 = np.array([KRAP,KSLC])
UA_4743 = np.array([KRAP,KDEN])
UA_5427 = np.array([KRAP,KORD])
AA_5765 = np.array([KRAP,KDFW])
QA___27 = np.array([YSSY,SCEL])
flights = ["DL_2599",
"DL_3710",
"UA_4743",
"UA_5427",
"AA_5765"]
fig, ax = plt.subplots(figsize=[11,8.5])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.NorthPolarStereo(central_longitude=-100.0)) # notice that
# here we are centering over CONUS!
ax.add_feature(feature = cfeature.COASTLINE)
ax.add_feature(feature = cfeature.STATES)
ax.plot(DL_2599[:,0], DL_2599[:,1], transform = ccrs.Geodetic())
ax.plot(DL_3710[:,0], DL_3710[:,1], transform = ccrs.Geodetic())
ax.plot(UA_4743[:,0], UA_4743[:,1], transform = ccrs.Geodetic())
ax.plot(UA_5427[:,0], UA_5427[:,1], transform = ccrs.Geodetic())
ax.plot(AA_5765[:,0], AA_5765[:,1], transform = ccrs.Geodetic())
ax.legend(flights)
ax.set_title("Polar Stereographic Projection")
ax.gridlines()
plt.show()
#
########################################################
/opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:825: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:877: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:944: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0:
########################################################
#
# Cartopy Example: Plotting A long flight route
#
# (0) Create our Flight Paths
#
# (1) Open a new default (single pannel) figure canvas
# (2) Reset the axis at (0,0) to be in Robinson Projection Coordinates
# (3) Set a title for the plot
# (4) Render continental boundaries.
# (5) Plot our flight paths using a Geodedic rendering!
# (6) Add Labels.
# (7) Add Grid Lines
#
# Bring out our flight paths!@
#
DL_2599 = np.array([KRAP,KMSP])
DL_3710 = np.array([KRAP,KSLC])
UA_4743 = np.array([KRAP,KDEN])
UA_5427 = np.array([KRAP,KORD])
AA_5765 = np.array([KRAP,KDFW])
QA___27 = np.array([YSSY,SCEL])
hirezSouthPS = ccrs.SouthPolarStereo()
hirezSouthPS._threshold = hirezSouthPS._threshold/100. #set finer mapping resolution
flights = ["QA___27"]
fig, ax = plt.subplots(figsize=[11,8.5])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=hirezSouthPS) # Deploy our custom-made projection!
ax.add_feature(feature = cfeature.COASTLINE)
ax.add_feature(feature = cfeature.BORDERS)
ax.plot(QA___27[:,0], QA___27[:,1], transform = ccrs.Geodetic())
ax.legend(flights)
ax.set_title("Polar Stereographic Projection")
ax.gridlines()
plt.show()
#
########################################################
/opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:825: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:877: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:944: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings.extend(multi_line_string) /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings.extend(multi_line_string)
And let's try another projection. The classic "Great Circle" routes are geodedic "as the crow flies" paths placed in a classic Mercator [conformal] projection.
This plot also uses cartopy.mpl.geoaxes.GeoAxes.set_extent() to fine tune the extents so give us a broader global perspective of our flights in the context of the big world.
########################################################
#
# Cartopy Example: Plotting A long flight route
#
# (0) Create our Flight Paths
#
# (1) Open a new default (single pannel) figure canvas
# (2) Reset the axis at (0,0) to be in Robinson Projection Coordinates
# (3) Set a title for the plot
# (4) Render continental boundaries.
# (5) Plot our flight paths using a Geodedic rendering!
# (6) Add Labels.
# (7) Fine Tune Map Extent
#
# Bring out our flight paths!@
#
DL_2599 = np.array([KRAP,KMSP])
DL_3710 = np.array([KRAP,KSLC])
UA_4743 = np.array([KRAP,KDEN])
UA_5427 = np.array([KRAP,KORD])
AA_5765 = np.array([KRAP,KDFW])
QA___27 = np.array([YSSY,SCEL])
hirezMercator = ccrs.Mercator(central_longitude = 180.0)
hirezMercator._threshold = hirezMercator._threshold/100. #set finer mapping resolution
flights = ["DL_2599",
"DL_3710",
"UA_4743",
"UA_5427",
"AA_5765",
"QA___27"]
fig, ax = plt.subplots(figsize=[11,8.5])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=hirezMercator) # Deploy our custom-made projection!
ax.add_feature(feature = cfeature.COASTLINE)
ax.add_feature(feature = cfeature.BORDERS)
ax.plot(DL_2599[:,0], DL_2599[:,1], transform = ccrs.Geodetic())
ax.plot(DL_3710[:,0], DL_3710[:,1], transform = ccrs.Geodetic())
ax.plot(UA_4743[:,0], UA_4743[:,1], transform = ccrs.Geodetic())
ax.plot(UA_5427[:,0], UA_5427[:,1], transform = ccrs.Geodetic())
ax.plot(AA_5765[:,0], AA_5765[:,1], transform = ccrs.Geodetic())
ax.plot(QA___27[:,0], QA___27[:,1], transform = ccrs.Geodetic())
ax.legend(flights)
ax.set_title("Mercator featuring 'Great Circle Route'")
# Fine tune the extents.
ax.set_extent([130, 360-50, -65, 50], # [min X, max X, min Y, max Y]
ccrs.PlateCarree()) # we need the Plate Carree
# lat/lon projextion into here
plt.show()
#
########################################################
/opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:825: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(multi_line_string) > 1: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:836: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings = list(multi_line_string) /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:836: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings = list(multi_line_string) /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:877: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. for line in multi_line_string: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:944: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. if len(p_mline) > 0: /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: Iteration over multi-part geometries is deprecated and will be removed in Shapely 2.0. Use the `geoms` property to access the constituent parts of a multi-part geometry. line_strings.extend(multi_line_string) /opt/miniconda3/lib/python3.9/site-packages/cartopy/crs.py:982: ShapelyDeprecationWarning: __len__ for multi-part geometries is deprecated and will be removed in Shapely 2.0. Check the length of the `geoms` property instead to get the number of parts of a multi-part geometry. line_strings.extend(multi_line_string)
If you need to break boundaries down in other countries to the equivalent of US States, one resource is the Natural Earth coverages. These provide a detailed set of natural and geopolitical boundaries.
Cartopy has an interface for it cartopy.feature.NaturalEarthFeature and it's demonstrated below for administrative units below the country level (e.g., US/Mexican/Australian States, French Departments, UK Counties, New Zealand Regions, etc).
Below that example is one that goes down one more level, down to US Counties (and only US Counties).
########################################################
#
# Cartopy Example: Plotting Domestic Flight Routes!
#
# (0) Fetch non-american provinces/states
#
# (1) Open a new default (single pannel) figure canvas (this time we're making it BIG!)
# (2) Reset the axis at (0,0) to be in Platte Carre Projection
# (3) Set a title for the plot
# (4) Pull the state-and-province level dataset from the Natural Earth Project as grey lines
# (5) Render our national borders and coastlines as black lines
fig, ax = plt.subplots(figsize=[17,11])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.PlateCarree())
ax.set_title("admin_1_states_provinces_lines")
# Create a feature for States/Admin 1 regions at 1:50m from Natural Earth
states_provinces = cfeature.NaturalEarthFeature(
category = 'cultural', # category physical or "cultural"
name = 'admin_1_states_provinces_lines', # name for the boudnary set
scale = '10m', # scale for the product (1:10m often works for me)
facecolor = 'none') # have a clear face color - keep it as 'none!'
ax.add_feature(states_provinces, edgecolor = 'darkgrey') # edgecolor will change the line colors
ax.add_feature(cfeature.BORDERS, edgecolor = 'black')
ax.add_feature(cfeature.COASTLINE, edgecolor = 'black')
plt.show()
#
########################################################
... and for US counties, including covering for large lake bodies...
########################################################
#
# Cartopy Example: Plotting Domestic Flight Routes!
#
# (0) Fetch non-american provinces/states
#
# (1) Open a new default (single pannel) figure canvas (this time we're making it BIG!)
# (2) Reset the axis at (0,0) to be in Platte Carre Projection
# (3) Set a title for the plot
# (4) Pull the state-and-province and county levels dataset from the Natural Earth Project as grey lines
# (5) Render various borders
fig, ax = plt.subplots(figsize=[17,11])
ax = plt.subplot(1, # First Row
1, # Second Column
1, # Start at Position 1 starting at one
projection=ccrs.PlateCarree())
ax.set_title("admin_2_counties")
# Create a feature for States/Admin 1 regions at 1:10m from Natural Earth
states_provinces = cfeature.NaturalEarthFeature(
category = 'cultural', # category physical or "cultural"
name = 'admin_1_states_provinces_lines', # name for the boudnary set
scale = '10m', # scale for the product (1:10m often works for me)
facecolor = 'none') # have a clear face color - keep it as 'none!'
# Create a feature for Cointies/Admin 2 regions at 1:10m from Natural Earth
#. the "_lake" will mask the lakes, id you need borders over lakes, remove it
us_counties = cfeature.NaturalEarthFeature(
category = 'cultural', # category physical or "cultural"
name = 'admin_2_counties_lakes', # name for the boudnary set
scale = '10m', # scale for the product (1:10m often works for me)
facecolor = 'none') # have a clear face color - keep it as 'none!'
ax.add_feature(us_counties,
edgecolor = 'grey') # edgecolor will change the line colors
ax.add_feature(states_provinces,
edgecolor = 'black',
linewidth = 1.5) # linethickness will change the line colors
ax.add_feature(cfeature.BORDERS,
edgecolor = 'black',
linewidth = 2)
ax.add_feature(cfeature.LAKES,
edgecolor = 'black',
linewidth = 1,
facecolor = "white") # facecolor="white" will mask lake areas
ax.add_feature(cfeature.COASTLINE,
edgecolor = 'black',
linewidth = 2)
ax.set_xlim((-93, -75))
ax.set_ylim(( 41, 50))
plt.show()
#
########################################################
There was once a very handy tool that would print version information of the Python version, operating system, and the versions of the Python packages you are running. It leveraged a "magic" command in IPython (which is what lies beneath Jupyter and JuptyerLab notebooks).
The developer has moved on to other things but the resource still has a narrow but ~militant~ ernest fan base (myself included). The original doesn't work with versions of Python above 3.7. So a couple of us wrote some patches to fix it. You can access my version below following these instructions.
!conda install -y -v git!pip install git+https://github.com/wjcapehart/version_informationFor people using Jupyter Lab, the interface does not play well with this "magic" function. It basically exports it as a JSON object with clicky expander buttons. However if you "Export" your notebook as an HTML or PDF file you should get a reasonable looking exported document.
################################################################
#
# Loading Version Information
#
%load_ext version_information
%version_information version_information, numpy, matplotlib, cartopy, shapely
#
################################################################
| Software | Version |
|---|---|
| Python | 3.9.7 64bit [Clang 11.1.0 ] |
| IPython | 7.30.1 |
| OS | macOS 12.0.1 x86\_64 i386 64bit |
| version_information | 1.0.3 |
| numpy | 1.19.5 |
| matplotlib | 3.5.0 |
| cartopy | 0.20.1 |
| shapely | 1.8.0 |
| Thu Dec 09 20:39:55 2021 MST | |