From 0332da45b1fec963ca9e9d1af39a68804b49e379 Mon Sep 17 00:00:00 2001 From: Jesse Pisel Date: Mon, 22 Jul 2019 21:53:11 -0600 Subject: [PATCH] Tutorial about how grid coordinates are generated (#192) Explain the intricacies of grid_coordinates, adjusting spacing/region, pixel registration, etc. Fixes #161 --- doc/index.rst | 1 + tutorials/grid_coordinates.py | 291 ++++++++++++++++++++++++++++++++++ 2 files changed, 292 insertions(+) create mode 100644 tutorials/grid_coordinates.py diff --git a/doc/index.rst b/doc/index.rst index c1fa75d43..981dc3402 100644 --- a/doc/index.rst +++ b/doc/index.rst @@ -29,6 +29,7 @@ :caption: User Guide sample_data/index.rst + tutorials/grid_coordinates.rst tutorials/trends.rst tutorials/decimation.rst tutorials/projections.rst diff --git a/tutorials/grid_coordinates.py b/tutorials/grid_coordinates.py new file mode 100644 index 000000000..cc5640190 --- /dev/null +++ b/tutorials/grid_coordinates.py @@ -0,0 +1,291 @@ +""" +.. _grid_coordinates: + +Grid Coordinates +==================== + +Creating the coordinates for regular grids in Verde is done using the +:func:`verde.grid_coordinates` function. It creates a set of regularly spaced points in +both the west-east and south-north directions, i.e. a two-dimensional spatial grid. These +points are then used by the Verde gridders to interpolate between data points. As such, all +`.grid` methods (like :meth:`verde.Spline.grid`) take as input the configuration parameters +for :func:`verde.grid_coordinates`. The grid can be specified either by the number of points +in each dimension (the shape) or by the grid node spacing. + + +""" + +import numpy as np +import matplotlib.pyplot as plt +import verde as vd +from matplotlib.patches import Rectangle + +######################################################################################## +# First let's create a region that is 1000 units west-east and 1000 units south-north, +# and we will set an initial spacing to 100 units. + +spacing = 100 +west, east, south, north = 0, 1000, 0, 1000 +region = (west, east, south, north) + +# create the grid coordinates +easting, northing = vd.grid_coordinates(region=region, spacing=spacing) + +######################################################################################## +# We can check the dimensions of the grid coordinates. The region is 1000 units and the +# spacing is 100 units, so the shape of the segments is 10x10. However, the number of +# grid nodes in this case is one more than the number of segments. So our grid +# coordinates have a shape of 11x11. +print(easting.shape, northing.shape) + +######################################################################################## +# Let's define two functions to visualize the region bounds and grid points + + +def plot_region(ax, region): + "Plot the region as a solid line." + west, east, south, north = region + ax.add_patch( + plt.Rectangle((west, south), east, north, fill=None, label="Region Bounds") + ) + + +def plot_grid(ax, coordinates, linestyles="dotted", region=None, pad=50, **kwargs): + "Plot the grid coordinates as dots and lines." + data_region = vd.get_region(coordinates) + ax.vlines( + coordinates[0][0], + ymin=data_region[2], + ymax=data_region[3], + linestyles=linestyles, + zorder=0, + ) + ax.hlines( + coordinates[1][:, 1], + xmin=data_region[0], + xmax=data_region[1], + linestyles=linestyles, + zorder=0, + ) + ax.scatter(*coordinates, **kwargs) + if pad: + padded = vd.pad_region(region, pad=pad) + plt.xlim(padded[:2]) + plt.ylim(padded[2:]) + + +######################################################################################## +# Visualize our region and grid coordinates using our functions + +plt.figure(figsize=(6, 6)) +ax = plt.subplot(111) +plot_region(ax=ax, region=region) +plot_grid( + ax=ax, + coordinates=(easting, northing), + region=region, + label="Square Region Grid Nodes", + marker=".", + color="black", + s=100, +) + +plt.xlabel("Easting") +plt.ylabel("Northing") +plt.legend(loc="upper center", bbox_to_anchor=(0.5, 1.15)) +plt.show() + +######################################################################################## +# Adjusting region boundaries when creating the grid coordinates +# -------------------------------------------------------------- +# +# Now let's change our spacing to 300 units. Because the range of the west-east and +# south-north boundaries are not multiples of 300, we must choose to change either: +# +# - the boundaries of the region in order to fit the spacing, or +# - the spacing in order to fit the region boundaries. +# +# We could tell :func:`verde.grid_coordinates` to adjust the region boundaries by +# passing ``adjust="region"``. + +spacing = 300 +region_easting, region_northing = vd.grid_coordinates( + region=region, spacing=spacing, adjust="region" +) +print(region_easting.shape, region_northing.shape) + +######################################################################################## +# With the spacing set at 300 units and a 4 by 4 grid of regular dimensions, +# :func:`verde.grid_coordinates` calculates the spatial location of each +# grid point and adjusts the region so that the maximum northing and maximum +# easting values are divisible by the spacing. In this example, both the easting and +# northing have 3 segments (4 nodes) that are each 300 units long, meaning the easting +# and northing span from 0 to 900. Both dimensions are divisible +# by 300. + +print(region_easting) +print(region_northing) + +###################################################################################### +# By default, if ``adjust`` is not assigned to ``"region"`` or ``"spacing"``, +# then :func:`verde.grid_coordinates` will adjust the spacing. With the adjust +# parameter set to ``spacing`` :func:`verde.grid_coordinates` creates grid nodes +# in a similar manner as when it adjusts the region. However, it doesn't readjust +# the region so that it is divisble by the spacing before creating the grid. +# This means the grid will have the same number of grid points no matter if +# the adjust parameter is set to ``region`` or ``spacing``. + +######################################################################################## +# Adjusting spacing when creating the grid +# ---------------------------------------- +# +# Now let's adjust the spacing of the grid points by passing ``adjust="spacing"`` +# to :func:`verde.grid_coordinates`. + +spacing_easting, spacing_northing = vd.grid_coordinates( + region=region, spacing=spacing, adjust="spacing" +) +print(spacing_easting.shape, spacing_northing.shape) + +###################################################################################### +# However the regular spacing between the grid points is no longer 300 units. + +print(spacing_easting) +print(spacing_northing) + +###################################################################################### +# Visualize the different adjustments +# ----------------------------------- +# +# Let's visualize the difference between the original region bounds, the +# adjusted region grid nodes, and the adjusted spacing grid nodes. Note the +# adjusted spacing grid nodes keep the original region, while the adjusted +# region grid nodes on the north and east side of the region have moved. + +plt.figure(figsize=(6, 6)) +ax = plt.subplot(111) +plot_region(ax=ax, region=region) +plot_grid( + ax=ax, + coordinates=(region_easting, region_northing), + region=region, + label="Adjusted Region Grid Nodes", + marker=">", + color="blue", + alpha=0.75, + s=100, +) +plot_grid( + ax=ax, + coordinates=(spacing_easting, spacing_northing), + region=region, + label="Adjusted Spacing Grid Nodes", + marker=">", + color="orange", + alpha=0.75, + s=100, +) +plt.xlabel("Easting") +plt.ylabel("Northing") +plt.legend(loc="upper center", bbox_to_anchor=(0.5, 1.18)) +plt.show() + +###################################################################################### +# Pixel Registration +# ------------------ +# +# Pixel registration locates the grid points in the middle of the grid segments +# rather than in the corner of each grid node. +# +# First, let's take our 1000x1000 region and use the 100 unit spacing from the first +# example and set the ``pixel_register`` parameter to ``True``. Without pixel +# registration our grid should have dimensions of 11x11. With pixel registration we +# expect the dimensions of the grid to be the dimensions of the non-registered grid +# minus one, or equal to the number of segments between the grid points in the +# non-registered grid (10x10). + +spacing = 100 +pixel_easting, pixel_northing = vd.grid_coordinates( + region=region, spacing=spacing, pixel_register=True +) +print(pixel_easting.shape, pixel_northing.shape) + +###################################################################################### +# And we can check the coordinates for the grid points with pixel registration. + +print(pixel_easting) +print(pixel_northing) + +###################################################################################### +# If we set ``pixel_register`` to ``False`` the function will return the grid +# coordinates of the nodes instead of pixel centers, resulting in an extra point in each direction. + +easting, northing = vd.grid_coordinates( + region=region, spacing=spacing, pixel_register=False +) +print(easting.shape, northing.shape) + +###################################################################################### +# Again we can check the coordinates for grid points with spacing adjustment. + +print(easting) +print(northing) + +###################################################################################### +# Lastly, we can visualize the pixel-registered grid points to see where they fall +# within the original region bounds. + +plt.figure(figsize=(6, 6)) +ax = plt.subplot(111) +plot_region(ax=ax, region=region) +plot_grid( + ax=ax, + coordinates=(pixel_easting, pixel_northing), + region=region, + label="Pixel Registered Grid Nodes", + marker=">", + color="blue", + alpha=0.75, + s=100, +) +plot_grid( + ax=ax, + coordinates=(easting, northing), + region=region, + label="Regular Registered Grid Nodes", + marker=">", + color="orange", + alpha=0.75, + s=100, +) + +plt.xlabel("Easting") +plt.ylabel("Northing") +plt.legend(loc="upper center", bbox_to_anchor=(0.5, 1.18)) +plt.show() + +###################################################################################### +# Extra Coordinates +# ----------------- +# +# In some cases, you might need an additional coordinate such as a height or a time +# that is associated with your coordinate grid. The ``extra_coords`` parameter +# in :func:`verde.grid_coordinates` creates an extra coordinate array that is the same +# shape as the coordinate grid, but contains a constant value. For example, let's +# add a constant height of 1000 units and time of 1 to our coordinate grid. + +easting, northing, height, time = vd.grid_coordinates( + region=region, spacing=spacing, extra_coords=[1000, 1] +) + +print(easting.shape, northing.shape, height.shape, time.shape) + +######################################################################################## +# And we can print the height array to verify that it is correct + +print(height) + +######################################################################################## +# And we can print the time array as well + +print(time)