172 lines
6.7 KiB
Plaintext
172 lines
6.7 KiB
Plaintext
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import cartopy.crs as ccrs\n",
|
|
"import cartopy.feature as cfeature\n",
|
|
"from cartopy.io.img_tiles import Stamen\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"from matplotlib.lines import Line2D as Line\n",
|
|
"from matplotlib.patheffects import Stroke\n",
|
|
"import numpy as np\n",
|
|
"import shapely.geometry as sgeom\n",
|
|
"from shapely.ops import transform as geom_transform\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The effect of badly referencing an ellipse\n",
|
|
"------------------------------------------\n",
|
|
"\n",
|
|
"This example demonstrates the effect of referencing your data to an incorrect\n",
|
|
"ellipse.\n",
|
|
"\n",
|
|
"First we define two coordinate systems - one using the World Geodetic System\n",
|
|
"established in 1984 and the other using a spherical globe. Next we extract\n",
|
|
"data from the Natural Earth land dataset and convert the Geodetic\n",
|
|
"coordinates (referenced in WGS84) into the respective coordinate systems\n",
|
|
"that we have defined. Finally, we plot these datasets onto a map assuming\n",
|
|
"that they are both referenced to the WGS84 ellipse and compare how the\n",
|
|
"coastlines are shifted as a result of referencing the incorrect ellipse."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"\n",
|
|
"def transform_fn_factory(target_crs, source_crs):\n",
|
|
" \"\"\"\n",
|
|
" Return a function which can be used by ``shapely.op.transform``\n",
|
|
" to transform the coordinate points of a geometry.\n",
|
|
"\n",
|
|
" The function explicitly *does not* do any interpolation or clever\n",
|
|
" transformation of the coordinate points, so there is no guarantee\n",
|
|
" that the resulting geometry would make any sense.\n",
|
|
"\n",
|
|
" \"\"\"\n",
|
|
" def transform_fn(x, y, z=None):\n",
|
|
" new_coords = target_crs.transform_points(source_crs,\n",
|
|
" np.asanyarray(x),\n",
|
|
" np.asanyarray(y))\n",
|
|
" return new_coords[:, 0], new_coords[:, 1], new_coords[:, 2]\n",
|
|
"\n",
|
|
" return transform_fn\n",
|
|
"\n",
|
|
"\n",
|
|
"def main():\n",
|
|
" # Define the two coordinate systems with different ellipses.\n",
|
|
" wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',\n",
|
|
" ellipse='WGS84'))\n",
|
|
" sphere = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',\n",
|
|
" ellipse='sphere'))\n",
|
|
"\n",
|
|
" # Define the coordinate system of the data we have from Natural Earth and\n",
|
|
" # acquire the 1:10m physical coastline shapefile.\n",
|
|
" geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))\n",
|
|
" dataset = cfeature.NaturalEarthFeature(category='physical',\n",
|
|
" name='coastline',\n",
|
|
" scale='10m')\n",
|
|
"\n",
|
|
" # Create a Stamen map tiler instance, and use its CRS for the GeoAxes.\n",
|
|
" tiler = Stamen('terrain-background')\n",
|
|
" fig = plt.figure()\n",
|
|
" ax = fig.add_subplot(1, 1, 1, projection=tiler.crs)\n",
|
|
" ax.set_title('The effect of incorrectly referencing the Solomon Islands')\n",
|
|
"\n",
|
|
" # Pick the area of interest. In our case, roughly the Solomon Islands, and\n",
|
|
" # get hold of the coastlines for that area.\n",
|
|
" extent = [155, 163, -11.5, -6]\n",
|
|
" ax.set_extent(extent, geodetic)\n",
|
|
" geoms = list(dataset.intersecting_geometries(extent))\n",
|
|
"\n",
|
|
" # Add the Stamen aerial imagery at zoom level 7.\n",
|
|
" ax.add_image(tiler, 7)\n",
|
|
"\n",
|
|
" # Transform the geodetic coordinates of the coastlines into the two\n",
|
|
" # projections of differing ellipses.\n",
|
|
" wgs84_geoms = [geom_transform(transform_fn_factory(wgs84, geodetic),\n",
|
|
" geom) for geom in geoms]\n",
|
|
" sphere_geoms = [geom_transform(transform_fn_factory(sphere, geodetic),\n",
|
|
" geom) for geom in geoms]\n",
|
|
"\n",
|
|
" # Using these differently referenced geometries, assume that they are\n",
|
|
" # both referenced to WGS84.\n",
|
|
" ax.add_geometries(wgs84_geoms, wgs84, edgecolor='white', facecolor='none')\n",
|
|
" ax.add_geometries(sphere_geoms, wgs84, edgecolor='gray', facecolor='none')\n",
|
|
"\n",
|
|
" # Create a legend for the coastlines.\n",
|
|
" legend_artists = [Line([0], [0], color=color, linewidth=3)\n",
|
|
" for color in ('white', 'gray')]\n",
|
|
" legend_texts = ['Correct ellipse\\n(WGS84)', 'Incorrect ellipse\\n(sphere)']\n",
|
|
" legend = ax.legend(legend_artists, legend_texts, fancybox=True,\n",
|
|
" loc='lower left', framealpha=0.75)\n",
|
|
" legend.legendPatch.set_facecolor('wheat')\n",
|
|
"\n",
|
|
" # Create an inset GeoAxes showing the location of the Solomon Islands.\n",
|
|
" sub_ax = fig.add_axes([0.7, 0.625, 0.2, 0.2],\n",
|
|
" projection=ccrs.PlateCarree())\n",
|
|
" sub_ax.set_extent([110, 180, -50, 10], geodetic)\n",
|
|
"\n",
|
|
" # Make a nice border around the inset axes.\n",
|
|
" effect = Stroke(linewidth=4, foreground='wheat', alpha=0.5)\n",
|
|
" sub_ax.outline_patch.set_path_effects([effect])\n",
|
|
"\n",
|
|
" # Add the land, coastlines and the extent of the Solomon Islands.\n",
|
|
" sub_ax.add_feature(cfeature.LAND)\n",
|
|
" sub_ax.coastlines()\n",
|
|
" extent_box = sgeom.box(extent[0], extent[2], extent[1], extent[3])\n",
|
|
" sub_ax.add_geometries([extent_box], ccrs.PlateCarree(), facecolor='none',\n",
|
|
" edgecolor='blue', linewidth=2)\n",
|
|
"\n",
|
|
" plt.show()\n"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"main()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"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.8.5"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 2
|
|
}
|