Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add a Derivative metagridder to predict derivatives #245

Draft
wants to merge 12 commits into
base: main
Choose a base branch
from
Prev Previous commit
Next Next commit
Rename Gradient to Derivative
  • Loading branch information
leouieda committed May 16, 2020
commit 61c90427be62e839db1d602598ec276f0a30660a
1 change: 1 addition & 0 deletions doc/api/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ Composite Estimators

Chain
Vector
Derivative

Model Selection
---------------
Expand Down
6 changes: 3 additions & 3 deletions examples/gradient.py
Original file line number Diff line number Diff line change
@@ -1,5 +1,5 @@
"""
Gradient calculations
Derivative calculations
=====================


Expand All @@ -15,8 +15,8 @@
).scatter(size=700, random_state=0)

spline = vd.Spline().fit((data.easting, data.northing), data.scalars)
east_deriv = vd.Gradient(spline, step=10, direction=(1, 0)).grid(spacing=50)
north_deriv = vd.Gradient(spline, step=10, direction=(0, 1)).grid(spacing=50)
east_deriv = vd.Derivative(spline, step=10, direction=(1, 0)).grid(spacing=50)
north_deriv = vd.Derivative(spline, step=10, direction=(0, 1)).grid(spacing=50)

fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(9, 10), sharex=True)

Expand Down
63 changes: 24 additions & 39 deletions examples/divergent_curl.py → examples/vector_derivatives.py
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
"""
Divergence and curl of vector fields
====================================
Derivatives of vector fields
============================
"""
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
Expand All @@ -11,78 +11,63 @@

# Fetch the wind speed data from Texas.
data = vd.datasets.fetch_texas_wind()
print(data.head())

# Separate out some of the data into utility variables
coordinates = (data.longitude.values, data.latitude.values)
region = vd.get_region(coordinates)
# Use a Mercator projection because Spline is a Cartesian gridder
projection = pyproj.Proj(proj="merc", lat_ts=data.latitude.mean())

# Split the data into a training and testing set. We'll fit the gridder on the training
# set and use the testing set to evaluate how well the gridder is performing.
# Split the data into a training and testing set. We'll fit the gridder on the
# training set and use the testing set to evaluate how well the gridder is
# performing.
train, test = vd.train_test_split(
projection(*coordinates),
(data.wind_speed_east_knots, data.wind_speed_north_knots),
random_state=2,
coordinates=projection(*coordinates),
data=(data.wind_speed_east_knots, data.wind_speed_north_knots),
test_size=0.1,
random_state=1,
)

# Chain together a blocked mean to avoid aliasing, a polynomial trend (Spline usually
# requires de-trended data), and finally a Spline for each component. Notice that
# BlockReduce can work on multicomponent data without the use of Vector.
chain = vd.Chain(
[
("trend", vd.Vector([vd.Trend(degree=1) for i in range(2)])),
(
"spline",
vd.Vector([vd.Spline(damping=1e-10, mindist=500e3) for i in range(2)]),
),
]
)
print(chain)
estimator = vd.Vector([vd.Spline(damping=1e-6, mindist=500e3) for i in range(2)])

# Fit on the training data
chain.fit(*train)
estimator.fit(*train)
# And score on the testing data. The best possible score is 1, meaning a perfect
# prediction of the test data.
score = chain.score(*test)
print("Cross-validation R^2 score: {:.2f}".format(score))
score = estimator.score(*test)
print("Validation score (R²): {:.2f}".format(score))

# Interpolate the wind speed onto a regular geographic grid and mask the data
# that are outside of the convex hull of the data points.
dims = ["latitude", "longitude"]
grid_full = chain.grid(region, spacing=20 / 60, projection=projection, dims=dims)
grid = vd.convexhull_mask(coordinates, grid=grid_full, projection=projection)

grid = estimator.grid(region, spacing=20 / 60, projection=projection, dims=dims)
grid = vd.convexhull_mask(coordinates, grid=grid, projection=projection)

spacing = 10 / 60
step = 0.5 * spacing * 100e3

east_derivs = vd.Gradient(chain, step=step, direction=(1, 0)).grid(
east_derivs = vd.Derivative(estimator, step=step, direction=(1, 0)).grid(
region, spacing=spacing, projection=projection, dims=dims
)
north_derivs = vd.Gradient(chain, step=step, direction=(0, 1)).grid(
north_derivs = vd.Derivative(estimator, step=step, direction=(0, 1)).grid(
region, spacing=spacing, projection=projection, dims=dims
)

divergence = east_derivs.east_component + north_derivs.north_component
curl = north_derivs.east_component - east_derivs.north_component


# Make maps of the original and gridded wind speed
plt.figure(figsize=(6, 6))
ax = plt.axes(projection=ccrs.Mercator())

curl.plot(ax=ax, transform=ccrs.PlateCarree())
divergence.plot(ax=ax, transform=ccrs.PlateCarree())

tmp = ax.quiver(
grid.longitude.values,
grid.latitude.values,
grid.east_component.values,
grid.north_component.values,
width=0.0015,
data.longitude.values,
data.latitude.values,
data.wind_speed_east_knots.values,
data.wind_speed_north_knots.values,
width=0.003,
scale=100,
color="tab:blue",
color="black",
transform=ccrs.PlateCarree(),
)

Expand Down
2 changes: 1 addition & 1 deletion verde/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,7 @@
pad_region,
longitude_continuity,
)
from .gradient import Gradient
from .derivative import Derivative
from .mask import distance_mask, convexhull_mask
from .utils import variance_to_weights, maxabs, grid_to_table
from .io import load_surfer
Expand Down
2 changes: 1 addition & 1 deletion verde/gradient.py → verde/derivative.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@
from .base.utils import check_data


class Gradient(BaseGridder):
class Derivative(BaseGridder):
"""
"""

Expand Down
32 changes: 16 additions & 16 deletions verde/tests/test_gradient.py → verde/tests/test_derivative.py
Original file line number Diff line number Diff line change
@@ -1,12 +1,12 @@
# pylint: disable=redefined-outer-name
"""
Test the Gradient estimator
Test the Derivative estimator
"""
import numpy as np
import numpy.testing as npt
import pytest

from ..gradient import Gradient
from ..derivative import Derivative
from ..trend import Trend
from ..vector import Vector
from ..chain import Chain
Expand Down Expand Up @@ -35,8 +35,8 @@ def test_gradient(polynomial):
coordinates, data, true_east_deriv, true_north_deriv = polynomial

trend = Trend(degree=2).fit(coordinates, data)
east_deriv = Gradient(trend, step=10, direction=(1, 0)).predict(coordinates)
north_deriv = Gradient(trend, step=10, direction=(0, 1)).predict(coordinates)
east_deriv = Derivative(trend, step=10, direction=(1, 0)).predict(coordinates)
north_deriv = Derivative(trend, step=10, direction=(0, 1)).predict(coordinates)

npt.assert_allclose(true_east_deriv, east_deriv, atol=1e-2)
npt.assert_allclose(true_north_deriv, north_deriv, atol=1e-2)
Expand All @@ -48,17 +48,17 @@ def test_gradient_fails_wrong_dimensions(polynomial):

trend = Trend(degree=2).fit(coordinates, data)
with pytest.raises(ValueError):
Gradient(trend, step=10, direction=(1, 0, 1)).predict(coordinates)
Derivative(trend, step=10, direction=(1, 0, 1)).predict(coordinates)
with pytest.raises(ValueError):
Gradient(trend, step=10, direction=(1, 0)).predict((coordinates[0],))
Derivative(trend, step=10, direction=(1, 0)).predict((coordinates[0],))


def test_gradient_grid(polynomial):
"Make sure the grid method works as expected"
coordinates, data, true_east_deriv = polynomial[:3]

trend = Trend(degree=2).fit(coordinates, data)
deriv = Gradient(trend, step=10, direction=(1, 0)).grid(spacing=50)
deriv = Derivative(trend, step=10, direction=(1, 0)).grid(spacing=50)

npt.assert_allclose(true_east_deriv, deriv.scalars.values, atol=1e-2)

Expand All @@ -70,22 +70,22 @@ def test_gradient_direction(polynomial):
for azimuth in np.radians(np.linspace(0, 360, 60)):
direction = (np.sin(azimuth), np.cos(azimuth))
true_deriv = true_east_deriv * direction[0] + true_north_deriv * direction[1]
deriv = Gradient(trend, step=10, direction=direction).predict(coordinates)
deriv = Derivative(trend, step=10, direction=direction).predict(coordinates)
npt.assert_allclose(true_deriv, deriv, atol=1e-2)


def test_gradient_fit(polynomial):
"Check that calling fit on the Gradient works as expected"
"Check that calling fit on the Derivative works as expected"
coordinates, data, true_east_deriv, true_north_deriv = polynomial

trend = Trend(degree=2)
east_deriv = (
Gradient(trend, step=10, direction=(1, 0))
Derivative(trend, step=10, direction=(1, 0))
.fit(coordinates, data)
.predict(coordinates)
)
north_deriv = (
Gradient(trend, step=10, direction=(0, 1))
Derivative(trend, step=10, direction=(0, 1))
.fit(coordinates, data)
.predict(coordinates)
)
Expand All @@ -95,14 +95,14 @@ def test_gradient_fit(polynomial):


def test_gradient_vector(polynomial):
"Make sure putting Gradients in a Vector works"
"Make sure putting Derivatives in a Vector works"
coordinates, data, true_east_deriv, true_north_deriv = polynomial

trend = Trend(degree=2)
gradient = Vector(
[
Gradient(trend, step=10, direction=(1, 0)),
Gradient(trend, step=10, direction=(0, 1)),
Derivative(trend, step=10, direction=(1, 0)),
Derivative(trend, step=10, direction=(0, 1)),
]
)
# This is wasteful because it fits the same trend twice so should not
Expand All @@ -115,11 +115,11 @@ def test_gradient_vector(polynomial):


def test_gradient_chain(polynomial):
"Make sure putting Gradients in a Chain works"
"Make sure putting Derivatives in a Chain works"
coordinates, data, true_east_deriv = polynomial[:3]

trend = Trend(degree=2)
gradient = Chain([("east", Gradient(trend, step=10, direction=(1, 0)))])
gradient = Chain([("east", Derivative(trend, step=10, direction=(1, 0)))])
gradient.fit(coordinates, data)
deriv = gradient.predict(coordinates)
npt.assert_allclose(true_east_deriv, deriv, atol=1e-2)