.ipynb .pdf

repository open issue suggest edit

Binder

Contents

Integration with the scientific Python ecosystem 🐍

In this tutorial, we’ll try out the integration between PyGMT and other common packages in the scientific Python ecosystem.

Besides pygmt, we’ll also be using:

• GeoPandas for managing geospatial tabular data

• Panel for interactive visualizations

• Xarray for managing n-dimensional labelled arrays

Plotting geospatial vector data with GeoPandas and PyGMT

We’ll extend the GeoPandas Mapping and Plotting Tools Examples to show how to create choropleth maps using PyGMT.

References:

• GeoPandas User Guide - https://geopandas.org/en/stable/docs/user_guide/

import pygmt
import geopandas as gpd

We’ll load sample data provided through the GeoPandas package and inspect the GeoDataFrame.

world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
world.head()

	pop_est	continent	name	iso_a3	gdp_md_est	geometry
0	920938	Oceania	Fiji	FJI	8374.0	MULTIPOLYGON (((180.00000 -16.06713, 180.00000...
1	53950935	Africa	Tanzania	TZA	150600.0	POLYGON ((33.90371 -0.95000, 34.07262 -1.05982...
2	603253	Africa	W. Sahara	ESH	906.5	POLYGON ((-8.66559 27.65643, -8.66512 27.58948...
3	35623680	North America	Canada	CAN	1674000.0	MULTIPOLYGON (((-122.84000 49.00000, -122.9742...
4	326625791	North America	United States of America	USA	18560000.0	MULTIPOLYGON (((-122.84000 49.00000, -120.0000...

Following the GeoPandas example, we’ll create a Choropleth map showing world population estimates, but will use PyGMT to plot the data using the Hammer projection.

# Calculate the populations in millions per capita
world = world[(world.pop_est>0) & (world.name!="Antarctica")]
world['pop_est'] = world.pop_est * 1e-6

# Find the range of data values for creating a colormap
cmap_bounds = pygmt.info(data=world['pop_est'], per_column=True)
cmap_bounds

array([1.40000000e-04, 1.37930277e+03])

Now, we’ll plot the data on a PyGMT figure, by creating a figure instance, laying down a basemap, plotting the GeoDataFrame, and adding a colorbar!

# Create an instance of the pygmt.Figure class
fig = pygmt.Figure()
# Create a colormap for the figure
pygmt.makecpt(cmap="bilbao", series=cmap_bounds)
# Create a basemap
fig.basemap(region="d", projection="H15c", frame=True)
# Plot the GeoDataFrame
# - Use `close=True` to specify that the polygons should be forced closed
# - Plot the polygon outlines with a 1 point, black pen
# - Set that the color should be based on the `pop_est` using the `color, `cmap`, and `aspatial` parameters
fig.plot(data=world, pen="1p,black", close=True, color="+z", cmap=True, aspatial="Z=pop_est")
# Add a colorbar
fig.colorbar(position="JMR", frame='a200+lPopulation (millions)')
# Display the output
fig.show()

/usr/share/miniconda3/envs/egu22pygmt/lib/python3.9/site-packages/geopandas/io/file.py:362: FutureWarning: pandas.Int64Index is deprecated and will be removed from pandas in a future version. Use pandas.Index with the appropriate dtype instead.
  pd.Int64Index,
plot [ERROR]: OGR parsing incomplete (is file missing OGR statements?) - abort

Interactive data visualization with Xarray, Panel, and PyGMT

In this section, we’ll create some interactive visualizations of oceanographic data!

We’ll use Panel, which is a Python library for connecting interactive widgets with plots! We’ll use Panel with PyGMT and xarray to visualize the objectively interpolated mean field for in-situ temperature from the World Ocean Atlas.

References:

• Temperature visualization based on https://rabernat.github.io/intro_to_physical_oceanography/02-c_ocean_temperature_salinity_stratification.html

• Interactive setup based on https://github.com/weiji14/30DayMapChallenge2021/blob/main/day25_interactive.py

• Data from the NOAA World Ocean Atlas, stored on the IRI Data Library at http://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NODC/.WOA09/.

import panel as pn
import xarray as xr
import pygmt
pn.extension()

# Download the dataset from the IRI Data Library
url = 'https://iridl.ldeo.columbia.edu/SOURCES/.NOAA/.NODC/.WOA09/.Grid-1x1/.Annual/.temperature/.t_an/data.nc'
netcdf_file = pygmt.which(fname=url, download=True)
woa_temp = xr.open_dataset(netcdf_file).isel(time=0)
woa_temp

<xarray.Dataset>
Dimensions:  (depth: 33, lat: 180, lon: 360)
Coordinates:
  * depth    (depth) float32 0.0 10.0 20.0 30.0 ... 4e+03 4.5e+03 5e+03 5.5e+03
  * lat      (lat) float32 -89.5 -88.5 -87.5 -86.5 -85.5 ... 86.5 87.5 88.5 89.5
  * lon      (lon) float32 0.5 1.5 2.5 3.5 4.5 ... 355.5 356.5 357.5 358.5 359.5
    time     datetime64[ns] 2008-01-01
Data variables:
    t_an     (depth, lat, lon) float32 ...

xarray.Dataset

• Dimensions:
  • depth: 33
  • lat: 180
  • lon: 360
• Coordinates: (4)
  • depth
    (depth)
    float32
    0.0 10.0 20.0 ... 5e+03 5.5e+03

    axis :

    Z

    positive :

    down

    description :

    Standard Depth Levels

    standard_name :

    depth

    long_name :

    depth

    gridtype :

    0

    units :

    m

    array([   0.,   10.,   20.,   30.,   50.,   75.,  100.,  125.,  150.,  200.,
            250.,  300.,  400.,  500.,  600.,  700.,  800.,  900., 1000., 1100.,
           1200., 1300., 1400., 1500., 1750., 2000., 2500., 3000., 3500., 4000.,
           4500., 5000., 5500.], dtype=float32)

  • lat
    (lat)
    float32
    -89.5 -88.5 -87.5 ... 88.5 89.5

    long_name :

    latitude

    bounds :

    lat_bnds

    pointwidth :

    1.0

    description :

    Latitude Range: -90 to 90 degrees with 0 at the Equator and increasing along the north.

    axis :

    Y

    standard_name :

    latitude

    gridtype :

    0

    units :

    degree_north

    array([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5, -80.5,
           -79.5, -78.5, -77.5, -76.5, -75.5, -74.5, -73.5, -72.5, -71.5, -70.5,
           -69.5, -68.5, -67.5, -66.5, -65.5, -64.5, -63.5, -62.5, -61.5, -60.5,
           -59.5, -58.5, -57.5, -56.5, -55.5, -54.5, -53.5, -52.5, -51.5, -50.5,
           -49.5, -48.5, -47.5, -46.5, -45.5, -44.5, -43.5, -42.5, -41.5, -40.5,
           -39.5, -38.5, -37.5, -36.5, -35.5, -34.5, -33.5, -32.5, -31.5, -30.5,
           -29.5, -28.5, -27.5, -26.5, -25.5, -24.5, -23.5, -22.5, -21.5, -20.5,
           -19.5, -18.5, -17.5, -16.5, -15.5, -14.5, -13.5, -12.5, -11.5, -10.5,
            -9.5,  -8.5,  -7.5,  -6.5,  -5.5,  -4.5,  -3.5,  -2.5,  -1.5,  -0.5,
             0.5,   1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,   9.5,
            10.5,  11.5,  12.5,  13.5,  14.5,  15.5,  16.5,  17.5,  18.5,  19.5,
            20.5,  21.5,  22.5,  23.5,  24.5,  25.5,  26.5,  27.5,  28.5,  29.5,
            30.5,  31.5,  32.5,  33.5,  34.5,  35.5,  36.5,  37.5,  38.5,  39.5,
            40.5,  41.5,  42.5,  43.5,  44.5,  45.5,  46.5,  47.5,  48.5,  49.5,
            50.5,  51.5,  52.5,  53.5,  54.5,  55.5,  56.5,  57.5,  58.5,  59.5,
            60.5,  61.5,  62.5,  63.5,  64.5,  65.5,  66.5,  67.5,  68.5,  69.5,
            70.5,  71.5,  72.5,  73.5,  74.5,  75.5,  76.5,  77.5,  78.5,  79.5,
            80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,  89.5],
          dtype=float32)

  • lon
    (lon)
    float32
    0.5 1.5 2.5 ... 357.5 358.5 359.5

    long_name :

    longitude

    bounds :

    lon_bnds

    pointwidth :

    1.0

    modulus :

    360.0

    description :

    Longitude Range: 0 to 360 degrees with 0 at Greenwich Meridian and increasing along the east.

    axis :

    X

    standard_name :

    longitude

    gridtype :

    1

    units :

    degree_east

    array([  0.5,   1.5,   2.5, ..., 357.5, 358.5, 359.5], dtype=float32)

  • time
    ()
    datetime64[ns]
    2008-01-01

    pointwidth :

    0

    gridtype :

    0

    array('2008-01-01T00:00:00.000000000', dtype='datetime64[ns]')

• Data variables: (1)
  • t_an
    (depth, lat, lon)
    float32
    ...

    cell_methods :

    time:mean within years time:mean over years

    grid_mapping :

    crs

    standard_name :

    sea_water_temperature

    long_name :

    Sea Water Temperature

    units :

    degrees Celsius

    description :

    The objectively interpolated mean field for Sea Water Temperature for annual time period at each one-degree square at each standard depth.

    [2138400 values with dtype=float32]

• Indexes: (3)
  • depth
    PandasIndex

    PandasIndex(Float64Index([   0.0,   10.0,   20.0,   30.0,   50.0,   75.0,  100.0,  125.0,
                   150.0,  200.0,  250.0,  300.0,  400.0,  500.0,  600.0,  700.0,
                   800.0,  900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0,
                  1750.0, 2000.0, 2500.0, 3000.0, 3500.0, 4000.0, 4500.0, 5000.0,
                  5500.0],
                 dtype='float64', name='depth'))

  • lat
    PandasIndex

    PandasIndex(Float64Index([-89.5, -88.5, -87.5, -86.5, -85.5, -84.5, -83.5, -82.5, -81.5,
                  -80.5,
                  ...
                   80.5,  81.5,  82.5,  83.5,  84.5,  85.5,  86.5,  87.5,  88.5,
                   89.5],
                 dtype='float64', name='lat', length=180))

  • lon
    PandasIndex

    PandasIndex(Float64Index([  0.5,   1.5,   2.5,   3.5,   4.5,   5.5,   6.5,   7.5,   8.5,
                    9.5,
                  ...
                  350.5, 351.5, 352.5, 353.5, 354.5, 355.5, 356.5, 357.5, 358.5,
                  359.5],
                 dtype='float64', name='lon', length=360))

• Attributes: (0)

# Make a static plot of sea surface temperature
fig = pygmt.Figure()
fig.grdimage(grid=woa_temp.t_an.sel(depth=0), cmap="vik", projection="R15c", frame=True)
fig.show()

grdinfo [WARNING]: Detected a data cube (/home/runner/work/egu22pygmt/egu22pygmt/book/data.nc) but -Q not set - skipping

# Make a panel widget for controlling the depth plotted
depth_slider = pn.widgets.DiscreteSlider(name='Depth (m)', options=woa_temp.depth.values.astype(int).tolist(), value=0)

# Make a function for plotting the depth slice with PyGMT

@pn.depends(depth=depth_slider)
def view(depth: int):
    fig = pygmt.Figure()
    pygmt.makecpt(cmap="vik", series=[-2,30])
    fig.grdimage(grid=woa_temp.t_an.sel(depth=depth), cmap=True, projection="R15c", frame=True)
    fig.colorbar(frame="a5")
    return fig

Make the interactive dashboard!

Now to put everything together! The ‘dashboard’ will be very simple. The ‘depth’ slider is placed next to the map using panel.Column. Selecting different depths will update the data plotted! Find out more at https://panel.holoviz.org/getting_started/index.html#using-panel.

Note: This is meant to run in a Jupyter lab/notebook environment. The grdinfo warning can be ignored.

pn.Column(depth_slider, view)

grdinfo [WARNING]: Detected a data cube (/home/runner/work/egu22pygmt/egu22pygmt/book/data.nc) but -Q not set - skipping

previous

Anatomy of a PyGMT figure

next

Making some Mars maps with pygmt