{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\nUsing weights in blocked means\n==============================\n\n:class:`verde.BlockReduce` is not able to output weights because we need to make\nassumptions about the reduction operation to know how to propagate uncertainties or\ncalculated weighted variances. That's why verde provides specialized reductions like\n:class:`verde.BlockMean`, which can calculate weights from input data in three ways:\n\n1. Using the variance of the data. This is the only possible option when no weights are\n provided.\n2. Using the uncertainty of the weighted mean propagated from the uncertainties in the\n data. In this case, we assume that the weights are ``1/uncertainty**2``.\n3. Using the weighted variance of the data.\n\nUsing the propagated uncertainties may be more adequate if your data is smooth in each\nblock (low variance) but have very different uncertainties. The propagation preserves a\nlow weight for data that have large uncertainties but don't vary much inside the block.\n\nThe weighted variance should be used when the data vary a lot in each block (high\nvariance) but have very similar uncertainties. The variance will be large when there is\na lot of variability in the data that isn't due to the uncertainties. This is also the\nbest choice if your data weights aren't ``1/uncertainty**2``.\n\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\nfrom matplotlib.colors import PowerNorm, LogNorm\nimport cartopy.crs as ccrs\nimport numpy as np\nimport verde as vd\n\n# We'll test this on the California vertical GPS velocity data because it comes with the\n# uncertainties\ndata = vd.datasets.fetch_california_gps()\ncoordinates = (data.longitude, data.latitude)\n\n# We'll calculate the mean on large blocks to show the effect of the different weighting\n# schemes\nspacing = 30 / 60\n# It's important that the weights are given as 1/sigma**2 for the uncertainty\n# propagation. In this case, you should not use verde.variance_to_weights because it\n# would normalize the weights.\nweights = 1 / data.std_up ** 2\nreducer = vd.BlockMean(spacing, center_coordinates=True)\n# First produce the weighted variance weights\nvariance_weights = reducer.filter(coordinates, data.velocity_up, weights)[-1]\n# And now produce the uncertainty propagation weights\nreducer.set_params(uncertainty=True)\nblock_coords, velocity, uncertainty_weights = reducer.filter(\n coordinates, data.velocity_up, weights\n)\n\n# Now we can plot the different weights side by side on Mercator maps\nfig, axes = plt.subplots(\n 1, 3, figsize=(13.5, 7), subplot_kw=dict(projection=ccrs.Mercator())\n)\ncrs = ccrs.PlateCarree()\ntitles = [\"Variance weights\", \"Uncertainty weights\"]\nweight_estimates = [variance_weights, uncertainty_weights]\nfor ax, title, w in zip(axes[:2], titles, weight_estimates):\n ax.set_title(title)\n # Plot the original data locations\n ax.plot(*coordinates, \".k\", transform=crs, markersize=0.5)\n # Plot the weights using a logarithmic color scale to highlight the\n # differences\n pc = ax.scatter(\n *block_coords, c=w, s=70, cmap=\"magma\", transform=crs, norm=LogNorm()\n )\n plt.colorbar(pc, ax=ax, orientation=\"horizontal\", pad=0.05)\n vd.datasets.setup_california_gps_map(ax)\n# Plot the original data uncertainties\nax = axes[2]\nax.set_title(\"Data uncertainty\")\n# Use a power law for the color scale because there is a lot of variability\n# Convert m/year to mm/year to have smaller values on the color bar\npc = ax.scatter(\n *coordinates,\n c=data.std_up * 1000,\n s=10,\n transform=crs,\n alpha=1,\n cmap=\"magma\",\n norm=PowerNorm(gamma=1 / 2)\n)\nplt.colorbar(pc, ax=ax, orientation=\"horizontal\", pad=0.05).set_label(\"mm/yr\")\nvd.datasets.setup_california_gps_map(ax)\nplt.tight_layout()\nplt.show()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.6"
}
},
"nbformat": 4,
"nbformat_minor": 0
}