Skip to content

Instantly share code, notes, and snippets.

@Guymer
Last active August 17, 2020 18:42
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save Guymer/004f602ea447fed829f50dbc6cb48f95 to your computer and use it in GitHub Desktop.
Save Guymer/004f602ea447fed829f50dbc6cb48f95 to your computer and use it in GitHub Desktop.
An example of some basic Cartopy functionality, that is optionally augmented by PyGuymer3 functionality.
#!/usr/bin/env python3
"""
This example scripts shows you how to do a few simple things using Cartopy. It
assumes a basic knowledge of MatPlotLib. This script has duplicated features too
as it also shows you how I personally choose to augment the Cartopy functions
with my own functions. This is optional and controlled by the "bellsWhistles"
Boolean variable. Finally, there are two functions in this script that contain
the "resolution" keyword which is set to "10m". If you want a cruder, more
jagged, plot that is generated quicker then you should set this to "110m". There
are only three possible choices: "10m"; "50m"; and "110m".
"""
# Import standard modules ...
import os
# Import special modules ...
try:
import cartopy
import cartopy.crs
except:
raise Exception("run \"pip install --user cartopy\"")
try:
import matplotlib
import matplotlib.pyplot
except:
raise Exception("run \"pip install --user matplotlib\"")
# Decide if you want bells and whistles ...
bellsWhistles = False
# Apply bells and whistles?
if bellsWhistles:
# Import my module ...
try:
import pyguymer3
except:
raise Exception("you need to have the Python module from https://github.com/Guymer/PyGuymer3 located somewhere in your $PYTHONPATH")
# Create figure, make the axes be the Robinson projection and make the extent be
# the whole world ...
# NOTE: Cartopy does the transformations (due to different map projections)
# automatically. All you need to tell it is the desired map projection of
# the output image (which is done here during the "ax = " call) and then
# you need to tell it the map projection of every piece of data that you
# draw. You can mix and match to your heart's content. You can even get it
# to connect points using either straight lines or great circles.
# * https://scitools.org.uk/cartopy/docs/latest/crs/projections.html -- supported projections
# * https://github.com/Guymer/fmc -- is an example that uses lon/lat coordinates
# * https://github.com/Guymer/hffl -- is an example that uses OSGB coordinates
fg = matplotlib.pyplot.figure(figsize = (9, 6), dpi = 300)
ax = matplotlib.pyplot.axes(projection = cartopy.crs.Robinson())
ax.set_global()
# Apply bells and whistles?
if bellsWhistles:
# Tell "pyguymer3" where to look for non-default background images ...
# NOTE: This folder is where you have placed the JSON and PNGs after
# following the tutorial on:
# * https://thomasguymer.co.uk/blog/2018/2018-01-15/
os.environ["CARTOPY_USER_BACKGROUNDS"] = "/path/to/NaturalEarthBackgroundImages"
# Add high-resolution background image just so that people can get their
# bearings ...
pyguymer3.add_map_background(ax, resolution = "large4096px")
else:
# Add low-resolution background image just so that people can get their
# bearings ...
ax.background_img()
# Draw the coastlines ...
ax.coastlines(resolution = "10m", color = "black", linewidth = 0.1)
# Draw the Tropic Of Cancer ...
# NOTE: See https://en.wikipedia.org/wiki/Tropic_of_Cancer
ax.plot(
[-180.0, 180.0],
[23.43662, 23.43662],
transform = cartopy.crs.PlateCarree(),
color = "black",
linewidth = 1.0,
linestyle = ":"
)
# Connect two points using a great circle ...
ax.plot(
[10.0, 80.0],
[30.0, 60.0],
transform = cartopy.crs.Geodetic(),
linewidth = 1.0,
color = "red"
)
# Connect the same two points using a straight-ish line ...
# NOTE: I don't know why you would want to do this - answers on a postcard
# please.
ax.plot(
[10.0, 80.0],
[30.0, 60.0],
transform = cartopy.crs.PlateCarree(),
linewidth = 1.0,
color = "green"
)
# Connect the same two points using a perfectly straight line ...
# NOTE: See https://stackoverflow.com/a/52861074
x1, y1 = cartopy.crs.Robinson().transform_point(10.0, 30.0, cartopy.crs.Geodetic())
x2, y2 = cartopy.crs.Robinson().transform_point(80.0, 60.0, cartopy.crs.Geodetic())
ax.plot(
[x1, x2],
[y1, y2],
transform = cartopy.crs.Robinson(),
linewidth = 1.0,
color = "blue"
)
# Find Shapefile containing all the country shapes ...
# NOTE: Cartopy automatically downloads free data from Natural Earth which I
# think is perfectly good enough for most applications, with the exception
# of Stansted being missing:
# * https://github.com/nvkelso/natural-earth-vector/issues/203
shape_file = cartopy.io.shapereader.natural_earth(resolution = "10m", category = "cultural", name = "admin_0_countries")
# Loop over records ...
for record in cartopy.io.shapereader.Reader(shape_file).records():
# Check if this country is Nepal ...
if record.attributes["NAME"] == "Nepal":
# Shade in Nepal ...
ax.add_geometries(
record.geometry,
cartopy.crs.PlateCarree(),
edgecolor = (1.0, 0.0, 0.0, 0.25),
facecolor = (1.0, 0.0, 0.0, 0.25),
linewidth = 0.5
)
# Save plot ...
fg.savefig("cartopy_example.png", bbox_inches = "tight", dpi = 300, pad_inches = 0.1)
# Apply bells and whistles?
if bellsWhistles:
# Shrink the PNG and remove EXIF data ...
pyguymer3.optimize_image("cartopy_example.png", strip = True)
# Clean up ...
matplotlib.pyplot.close("all")
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment