.. _sphx_glr_gallery_effects_of_the_ellipse.py:


The effect of badly referencing an ellipse
------------------------------------------

This example demonstrates the effect of referencing your data to an incorrect
ellipse.

First we define two coordinate systems - one using the World Geodetic System
established in 1984 and the other using a spherical globe. Next we extract
data from the Natural Earth land dataset and convert the Geodetic
coordinates (referenced in WGS84) into the respective coordinate systems
that we have defined. Finally, we plot these datasets onto a map assuming
that they are both referenced to the WGS84 ellipse and compare how the
coastlines are shifted as a result of referencing the incorrect ellipse.





.. image:: /gallery/images/sphx_glr_effects_of_the_ellipse_001.png
    :align: center





.. code-block:: python


    import cartopy.crs as ccrs
    import cartopy.feature as cfeature
    from cartopy.io.img_tiles import StamenTerrain
    import matplotlib.pyplot as plt
    from matplotlib.lines import Line2D as Line
    from matplotlib.patheffects import Stroke
    import numpy as np
    import shapely.geometry as sgeom
    from shapely.ops import transform as geom_transform


    def transform_fn_factory(target_crs, source_crs):
        """
        Return a function which can be used by ``shapely.op.transform``
        to transform the coordinate points of a geometry.

        The function explicitly *does not* do any interpolation or clever
        transformation of the coordinate points, so there is no guarantee
        that the resulting geometry would make any sense.

        """
        def transform_fn(x, y, z=None):
            new_coords = target_crs.transform_points(source_crs,
                                                     np.asanyarray(x),
                                                     np.asanyarray(y))
            return new_coords[:, 0], new_coords[:, 1], new_coords[:, 2]

        return transform_fn


    def main():
        # Define the two coordinate systems with different ellipses.
        wgs84 = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
                                                  ellipse='WGS84'))
        sphere = ccrs.PlateCarree(globe=ccrs.Globe(datum='WGS84',
                                                   ellipse='sphere'))

        # Define the coordinate system of the data we have from Natural Earth and
        # acquire the 1:10m physical coastline shapefile.
        geodetic = ccrs.Geodetic(globe=ccrs.Globe(datum='WGS84'))
        dataset = cfeature.NaturalEarthFeature(category='physical',
                                               name='coastline',
                                               scale='10m')

        # Create a Stamen map tiler instance, and use its CRS for the GeoAxes.
        tiler = StamenTerrain()
        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1, projection=tiler.crs)
        ax.set_title('The effect of incorrectly referencing the Solomon Islands')

        # Pick the area of interest. In our case, roughly the Solomon Islands, and
        # get hold of the coastlines for that area.
        extent = [155, 163, -11.5, -6]
        ax.set_extent(extent, geodetic)
        geoms = list(dataset.intersecting_geometries(extent))

        # Add the Stamen aerial imagery at zoom level 7.
        ax.add_image(tiler, 7)

        # Transform the geodetic coordinates of the coastlines into the two
        # projections of differing ellipses.
        wgs84_geoms = [geom_transform(transform_fn_factory(wgs84, geodetic),
                                      geom) for geom in geoms]
        sphere_geoms = [geom_transform(transform_fn_factory(sphere, geodetic),
                                       geom) for geom in geoms]

        # Using these differently referenced geometries, assume that they are
        # both referenced to WGS84.
        ax.add_geometries(wgs84_geoms, wgs84, edgecolor='white', facecolor='none')
        ax.add_geometries(sphere_geoms, wgs84, edgecolor='gray', facecolor='none')

        # Create a legend for the coastlines.
        legend_artists = [Line([0], [0], color=color, linewidth=3)
                          for color in ('white', 'gray')]
        legend_texts = ['Correct ellipse\n(WGS84)', 'Incorrect ellipse\n(sphere)']
        legend = ax.legend(legend_artists, legend_texts, fancybox=True,
                           loc='lower left', framealpha=0.75)
        legend.legendPatch.set_facecolor('wheat')

        # Create an inset GeoAxes showing the location of the Solomon Islands.
        sub_ax = fig.add_axes([0.7, 0.625, 0.2, 0.2],
                              projection=ccrs.PlateCarree())
        sub_ax.set_extent([110, 180, -50, 10], geodetic)

        # Make a nice border around the inset axes.
        effect = Stroke(linewidth=4, foreground='wheat', alpha=0.5)
        sub_ax.outline_patch.set_path_effects([effect])

        # Add the land, coastlines and the extent of the Solomon Islands.
        sub_ax.add_feature(cfeature.LAND)
        sub_ax.coastlines()
        extent_box = sgeom.box(extent[0], extent[2], extent[1], extent[3])
        sub_ax.add_geometries([extent_box], ccrs.PlateCarree(), facecolor='none',
                              edgecolor='blue', linewidth=2)

        plt.show()


    if __name__ == '__main__':
        main()

**Total running time of the script:** ( 0 minutes  0.683 seconds)



.. container:: sphx-glr-footer


  .. container:: sphx-glr-download

     :download:`Download Python source code: effects_of_the_ellipse.py <effects_of_the_ellipse.py>`



  .. container:: sphx-glr-download

     :download:`Download Jupyter notebook: effects_of_the_ellipse.ipynb <effects_of_the_ellipse.ipynb>`

.. rst-class:: sphx-glr-signature

    `Generated by Sphinx-Gallery <https://sphinx-gallery.readthedocs.io>`_