Examples of using Cartopy with Julia

Julia Cartopy Examples

Last Update: 2024-09-13

Installation
3D Earth observed from infinite distance
3D Earth observed from finite altitude
2D world map
2D world map with stock image
Draw points, lines, & great circles
Add features
Overlay images
Overlay vectors
Threshold smoothing

Cartopy Installation

In the Julia command line: using Conda; Conda.add("Cartopy")
In computer's command line (Terminal on Mac, Command Prompt on Windows): conda install -c conda-forge cartopy

Python → Julia PyPlot Translation

When using Python modules (of which Cartopy is an example) in Julia rather than in Python, the syntax generally changes as follows:

Python & Julia 1.0: ax.stock_img()
Julia ≤0.7: ax[:stock_img]()

3D Earth observed from infinite distance

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
feature = pyimport("cartopy.feature")

lat = 15 # °
lon = -50 # °
prj = ccrs.Orthographic(lon, lat)

ax = subplot(projection=prj)
ax.set_global() # make the map global rather than have it zoom in to the extents of any plotted data

ax.add_feature(feature.OCEAN, color="navy")
ax.add_feature(feature.LAND, color="lightgray")

axis("off")

3D Earth observed from finite altitude

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
feature = pyimport("cartopy.feature")

lat = 38 # °
lon = 14 # °
alt = 100e3 # meters

prj = ccrs.NearsidePerspective(central_latitude=lat, central_longitude=lon, satellite_height=alt)

ax = subplot(projection=prj)
ax.set_global()
ax.gridlines(lw=0.5, ls="--")

ax.add_feature(feature.OCEAN, color="navy")
ax.add_feature(feature.LAND, color="lightgray")

axis("off")

2D world map

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
feature = pyimport("cartopy.feature")

prj = ccrs.PlateCarree()

ax = subplot(projection=prj)
ax.set_global()

ax.add_feature(feature.OCEAN, color="navy")
ax.add_feature(feature.LAND, color="lightgray")

axis("off")

2D world map with stock image

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")

prj = ccrs.PlateCarree()

ax = subplot(projection=prj)
ax.set_global()
ax.stock_img()

axis("off")

Draw points, lines, & great circles

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")

prj = ccrs.Robinson()

ax = subplot(projection=prj)

ax.set_global()
ax.stock_img()
ax.coastlines()

tfm1 = ccrs.PlateCarree()
tfm2 = ccrs.Geodetic()

scatter(-0.1, 51.5, transform=tfm1, zorder=4, s=40, linewidth=1.5, edgecolor="k", color="yellow")
scatter(132, 43.2, transform=tfm1, zorder=4, s=80, linewidth=2, edgecolor="b", color="c")
scatter(-58.4, -34.6, transform=tfm1, zorder=4, s=60, linewidth=2, edgecolor="g", color="orange")

plot([-0.1, 132], [51.5, 43.2], transform=tfm1, linewidth=3, "r") # Straight line on 2D map
plot([-58.4, 132], [-34.6, 43.2], transform=tfm2, linewidth=3, "m") # Curved line (great circle) on 2D map

title("Global Map with Points, Lines, and Great Circles")

Add features

Such as state borders, rivers, coastlines, etc.

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
feature = pyimport("cartopy.feature")

prj = ccrs.PlateCarree()

ax = subplot(projection=prj)
ax.set_extent([80, 170, -45, 30])

ax.stock_img()

# Create a feature for States/Admin 1 regions at 1:50m from Natural Earth
states_provinces = feature.NaturalEarthFeature(
        category="cultural",
        name="admin_1_states_provinces_lines",
        scale="50m",
        facecolor="none"
)

ax.add_feature(feature.LAND)
ax.add_feature(feature.COASTLINE)
ax.add_feature(states_provinces, edgecolor="gray")

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
feature = pyimport("cartopy.feature")

prj = ccrs.PlateCarree()

ax = subplot(projection=prj)
ax.set_extent([-20, 60, -40, 40])

ax.add_feature(feature.LAND)
ax.add_feature(feature.OCEAN)
ax.add_feature(feature.COASTLINE)
ax.add_feature(feature.BORDERS, linestyle=":")
ax.add_feature(feature.LAKES, alpha=0.5)
ax.add_feature(feature.RIVERS)

Overlay images

using PyPlot, PyCall
cartopy = pyimport("cartopy")

fig = figure(figsize=(8, 12))

fname = joinpath(cartopy.config["repo_data_dir"], "raster", "sample", "Miriam.A2012270.2050.2km.jpg")
img_extent = (-120.6766, -106.3210, 13.2301, 30.7669)
img = imread(fname)

prj = cartopy.crs.PlateCarree()
ax = subplot(projection=prj)

title("Hurricane Miriam from the Aqua/MODIS satellite\n2012 09/26/2012 20:50 UTC")

xlim(img_extent[1]-1, img_extent[2]+1) # set a margin around the data

# add the image. Because this image was a tif, the "origin" of the image is in the upper left corner
tfm1 = cartopy.crs.PlateCarree()
imshow(img, origin="upper", extent=img_extent, transform=tfm1)
ax.coastlines(resolution="50m", color="k", linewidth=1)

# mark a known place to help us geo-locate ourselves
tfm2 = cartopy.crs.Geodetic()
plot(-117.1625, 32.715, "bo", markersize=7, transform=tfm2)
text(-117, 33, "San Diego", transform=ccrs.Geodetic())

Overlay vectors

using PyCall, PyPlot
ccrs = pyimport("cartopy.crs")
feature = pyimport("cartopy.feature")

"""
Returns `(x, y, u, v)` of some vector data
computed mathematically. 
"""
function sample_data(shape=(20, 30))
    x = range(311.9, stop=391.1, length=shape[2])
    y = range(-23.6, stop=24.8, length=shape[1])

    x2d, y2d = repeat(x', length(y), 1), repeat(y, 1, length(x))
    u = 10(2cos.(2deg2rad.(x2d) .+ 3deg2rad.(y2d .+ 30)).^2)
    v = 20cos.(6deg2rad.(x2d))

    return x, y, u, v
end

prj = ccrs.Orthographic(-10, 45)
ax = subplot(projection=prj)

ax.set_global()
ax.gridlines()
ax.add_feature(feature.OCEAN, zorder=0)
ax.add_feature(feature.LAND, zorder=0, edgecolor="k")

# The transform will be a rotated pole CRS, meaning that the vectors will be unevenly spaced in regular PlateCarree space.
tfm = ccrs.RotatedPole(pole_longitude=177.5, pole_latitude=37.5)

x, y, u, v = sample_data()
ax.quiver(collect(x), collect(y), u, v, transform=tfm)

Multiple Projections, Multiple Subplots

Smoothing great circle lines by adjusting projection threshold.

using PyPlot, PyCall
ccrs = pyimport("cartopy.crs")
feature = pyimport("cartopy.feature")

coord1 = (120, 45)
coord2 = (-58, -35)

fig = figure(figsize=(8,8))

tfm1 = ccrs.PlateCarree()
tfm2 = ccrs.Geodetic()

prj1 = ccrs.PlateCarree()
prj1.threshold = 2
ax1 = subplot(211, projection=prj1)

prj2 = ccrs.PlateCarree()
prj2.threshold = 1e-2
ax2 = subplot(212, projection=prj2)

for ax in [ax1, ax2]
    ax.set_global()

    ax.add_feature(feature.OCEAN, color="white")
    ax.add_feature(feature.LAND, color="lightgray")

    ax.scatter(coord1[1], coord1[2], transform=tfm1, zorder=4, color="k", s=40)
    ax.scatter(coord2[1], coord2[2], transform=tfm1, zorder=4, color="k", s=40)

    ax.plot([coord1[1], coord2[1]], [coord1[2], coord2[2]], transform=tfm2, linewidth=3, "k")
end

jpwspicer/Julia Cartopy Examples.md