Extract and plot data along a profile#
The module polartoolkit.profile provides convenient tools to extract data along a profile, and automatically create plots of the results. These plots can optionally include a map, a cross section of earth layers, and a profile of other datasets.
There are 3 methods to define the location of the profile: “Points”, “Shapefile”, and “Polyline”
[2]:
%%capture
%load_ext autoreload
%autoreload 2
from polartoolkit import fetch, profile, utils
Method 1: “Points”#
define a start and end point for a profile, in EPSG:3031 projection. See {doc}../tips for a map with EPSG gridlines. Sample bedmap2 layers between these points, and plot the cross section
[3]:
a = (-1200e3, -1400e3)
b = (1000e3, 1400e3)
[4]:
# this functions returns the figure as well as 2 dataframe with the sampled layers and
# sampled data
fig, df_layers, df_data = profile.plot_profile(
method="points",
start=a,
stop=b,
)
fig.show()
Method 2: “Shapefile”:#
Instead of using a line defined by 2 points, we can sample along the path of a shapefile. Provide a file name, or load a already included shapefile with the function fetch.
Here is a curved path across the Ross Ice Shelf from Mulock Glacier to the ice front through Discovery Deep.
[5]:
fig, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp("Disco_deep_transect"),
fillnans=True,
)
fig.show()
Method 3: “Polyline”:#
Alternatively, load an interactive map, click to create a line, and use that line to create the path.
Note: this requires a few optional dependencies. This can be install with: pip install polartoolkit[interactive] or installed individually with: mamba install geoviews cartopy ipyleaflet ipython
The below cell will display an interactive map of Antarctica. Drag or zoom to your region of interest and use the Draw a polyline button to create the profile you want.
[10]:
lines = profile.draw_lines()
[11]:
# extract the vertices of the line
df = utils.shapes_to_df(lines)
print(df)
lon lat shape_num x y
0 -57.508101 -74.098180 0 -1.466368e+06 933888.012717
1 -61.032339 -75.595416 0 -1.376256e+06 761856.010111
2 -67.463212 -78.042204 0 -1.204224e+06 499711.993138
3 -81.408775 -79.429845 0 -1.138688e+06 172031.988060
4 -98.315126 -80.117920 0 -1.064960e+06 -155648.000321
5 -119.010360 -80.537226 0 -9.011200e+05 -499711.970129
6 -136.347887 -80.952321 0 -6.799360e+05 -712703.983999
[12]:
# use the line to create a profile, and resample it to get 1000 points
fig, _, _ = profile.plot_profile(
method="polyline",
polyline=df,
num=1000,
)
fig.show()
Add a map#
show the profile location by adding a map with a default background of satellite imagery
[8]:
fig, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp("Disco_deep_transect"),
add_map=True,
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Add a graph of data#
Sample and plot additional datasets along the same profile with data_dict. By default this includes gravity and magnetics.
[9]:
fig, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp("Disco_deep_transect"),
add_map=True,
data_dict="default",
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
save=True, # save the plot to use in README
path="../cover_fig.png",
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Change profile length#
clip the profile (either end) based on distance
[10]:
fig, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp("Disco_deep_transect"),
data_dict="default",
add_map=True,
clip=True,
min_dist=50e3,
max_dist=200e3,
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Change sampling resolution#
increase the resolution of the sampling. If using a profile defined by 2 points, use parameter “num”, which defaults to 100 points along the line.
[11]:
fig1, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp(
"Roosevelt_Island",
),
num=30,
data_dict="default",
add_map=True,
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig1.show()
fig2, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp(
"Roosevelt_Island",
),
num=200,
data_dict="default",
add_map=True,
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig2.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
gmtset [WARNING]: Representation of font type not recognized. Using default.
Alter map properies#
Use map_background to change the grid shown on the map. Here, instead of the default imagery, we will show the surface topography from bedmap2.
Also, use subplot_orientation to stack the plots vertically
[12]:
fig, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp("Disco_deep_transect"),
data_dict="default",
add_map=True,
map_background=fetch.bedmap2("surface", fill_nans=True),
subplot_orientation="vertical",
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Use map_cmap to change the maps colorscale
[13]:
fig, _, _ = profile.plot_profile(
method="shapefile",
shapefile=fetch.sample_shp("Disco_deep_transect"),
data_dict="default",
add_map=True,
map_background=fetch.bedmap2("surface", fill_nans=True),
map_cmap="viridis",
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Use map_buffer to change the level of zoom on the map. This is a percentage of total line distance, which defaults to 0.3 (130%) of profile length. So if the profile is 100km long, the default map will be 130km along the longest axis.
[14]:
# define new profile endpoints
a = (-590e3, -1070e3)
b = (-100e3, -545e3)
fig1, _, _ = profile.plot_profile(
method="points",
start=a,
stop=b,
data_dict="default",
add_map=True,
map_buffer=0.1,
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig1.show()
fig2, _, _ = profile.plot_profile(
method="points",
start=a,
stop=b,
data_dict="default",
add_map=True,
map_buffer=0.8,
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig2.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
gmtset [WARNING]: Representation of font type not recognized. Using default.
Alter Cross-section properties#
Use layer_buffer to alter the amount of whitespace above and below the cross-section layers, this defaults to 0.1 (10% of data spread added above/below)
[15]:
fig, _, _ = profile.plot_profile(
method="points",
start=a,
stop=b,
data_dict="default",
add_map=True,
layer_buffer=1.0,
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Alter Cross-section properties#
Use data_buffer to alter the whitespace above and below data graph, defaults to 0.1 (10% of data spread added above/below)
[16]:
fig, _, _ = profile.plot_profile(
method="points",
start=a,
stop=b,
data_dict="default",
add_map=True,
data_buffer=0.8,
legend_loc="JTR+jTR+o.2c/-6c", # change the default legend location with GMT format
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Custom layers, datasets, and colors#
The input for both the data and the layers are nested dictionaries, where each dictionary entry contains keys: ‘name’, ‘grid’, and ‘color’.
Use the function, profile.make_data_dict() to help create these dictionaries. It takes 3 inputs, which are: * a list of names * a list of grids either geotiff/netcdf file names or xarray.DataArrays * a list of colors.
For example, if you have a netcdf file ‘ice_velocity.nc’ which you want to plot as orange, you can make that into a dictionary with: profile.make_data_dict(['ice velocity'], ['ice_velocity.nc'], ['orange'])
Optionally, you can use Pooch.fetch to help download and store dataset from urls, and easily load them.
See the functions in the module polartoolkit.fetch for examples and feel free to make a PR and add your datasets there.
[17]:
# note, these datafile can be plotted at higher resolution by adding parameter 'spacing'
# to the 'fetch' function:
# fetch.gravity(version="antgg-update", anomaly_type="BA", spacing=1e3) will use a 1km
# version of the Bouguer gravity anomaly.
data_dict = profile.make_data_dict(
["Bouguer gravity", "REMA surface elevation"],
[
fetch.gravity(
version="antgg-update",
anomaly_type="BA",
),
fetch.rema(version="1km"),
],
["purple", "red"],
)
# get default bedmap2 layers
layers_dict = profile.default_layers(version="bedmap2")
# add dictionary entry of extra layer 'basement'
layers_dict["basement"] = {}
layers_dict["basement"]["name"] = "basement"
layers_dict["basement"]["grid"] = fetch.basement()
layers_dict["basement"]["color"] = "chocolate"
a = (-590e3, -1070e3)
b = (-100e3, -545e3)
fig, _, _ = profile.plot_profile(
"points",
start=a,
stop=b,
add_map=True,
data_dict=data_dict,
layers_dict=layers_dict,
legend_loc="JTR+jTR+o.2c/-6.3c", # change the default legend location
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
Control wheather to fill gaps in the cross-section layers. Note that the additional layer ‘basement’, doesn’t extend past the groundingline. By default, NaN’s in any layer are set equal to the layer above, causing the vertical line at ~100km.
[18]:
fig, _, _ = profile.plot_profile(
"points",
start=a,
stop=b,
add_map=True,
data_dict=data_dict,
layers_dict=layers_dict,
fillnans=False,
legend_loc="JTR+jTR+o.2c/-6.3c", # change the default legend location
)
fig.show()
gmtset [WARNING]: Representation of font type not recognized. Using default.
[ ]: