Source code for epygram._plugins.with_cartopy.D3ProjectedGeometry

#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Copyright (c) Météo France (2014-)
# This software is governed by the CeCILL-C license under French law.
# http://www.cecill.info
"""
Extend D3Geometry with plotting methods using cartopy.
"""
from __future__ import print_function, absolute_import, unicode_literals, division

import numpy

from cartopy import crs as ccrs

import footprints

from epygram import epygramError

epylog = footprints.loggers.getLogger(__name__)


def activate():
    """Activate extension."""
    from . import __name__ as plugin_name
    from epygram._plugins.util import notify_doc_requires_plugin
    notify_doc_requires_plugin([default_cartopy_CRS, get_cartopy_extent],
                               plugin_name)
    from epygram.geometries.D3Geometry import D3ProjectedGeometry
    # defaults arguments for cartopy plots
    D3ProjectedGeometry.default_cartopy_CRS = default_cartopy_CRS
    D3ProjectedGeometry.get_cartopy_extent = get_cartopy_extent


def default_cartopy_CRS(self):
    """
    Create a cartopy.crs appropriate to the Geometry.
    """
    if self.name == 'lambert':
        if self.secant_projection:
            lat_0 = (self.projection['secant_lat1'].get('degrees') +
                     self.projection['secant_lat2'].get('degrees')) / 2.
            secant_lats = (self.projection['secant_lat1'].get('degrees'),
                           self.projection['secant_lat2'].get('degrees'))
        else:
            lat_0 = self.projection['reference_lat'].get('degrees')
            secant_lats = (lat_0, lat_0)
        crs = ccrs.LambertConformal(
            central_longitude=self.projection['reference_lon'].get('degrees'),
            central_latitude=lat_0,
            standard_parallels=secant_lats)
    elif self.name == 'mercator':
        if self.secant_projection:
            lat = 'secant_lat'
        else:
            lat = 'reference_lat'
        crs = ccrs.Mercator(
            central_longitude=self._center_lon.get('degrees'),
            latitude_true_scale=self.projection[lat].get('degrees'))
    elif self.name == 'polar_stereographic':
        if self.secant_projection:
            lat = 'secant_lat'
        else:
            lat = 'reference_lat'
        crs = ccrs.Stereographic(
            central_latitude=numpy.copysign(90., self.projection[lat].get('degrees')),
            central_longitude=self.projection['reference_lon'].get('degrees'),
            true_scale_latitude=self.projection[lat].get('degrees'))
    elif self.name == 'space_view':
        if self.projection['satellite_lat'].get('degrees') == 0.:
            crs = ccrs.Geostationary(
                central_longitude=self.projection['satellite_lon'].get('degrees'),
                satellite_height=self.projection['satellite_height'])
        else:
            raise NotImplementedError("Implementatation to be tested")
            crs = ccrs.NearsidePerspective(
                central_longitude=self.projection['satellite_lon'].get('degrees'),
                central_latitude=self.projection['satellite_lat'].get('degrees'),
                satellite_height=self.projection['satellite_height'])
    else:
        raise epygramError("Projection name unknown.")
    return crs


def get_cartopy_extent(self, subzone=None):
    """
    Gets the extension of the geometry in the default crs
    :param subzone: defines the LAM subzone to be included, in LAM case,
                    among: 'C', 'CI'.
    :return: (xmin, xmax, ymin, ymax) as expected by matplotlib's ax.set_extent
    """
    # Actually, this method transforms projected coordinates from pyproj to
    # projected coordinates from cartopy's CRS
    # Both make use of proj4 but the transformation is not exactly the same.
    # Easiest method is to transform pyproj coordinates into lat/lon then
    # into crs coordinates. But this method does not work when the point
    # does not have lat/lon. This is the case in space_view where corners
    # can be outside the sphere.
    crs = self.default_cartopy_CRS()
    if self.name != 'space_view':
        (llcrnrlon, llcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ll'])
        (urcrnrlon, urcrnrlat) = self.ij2ll(*self.gimme_corners_ij(subzone)['ur'])
        ll = crs.transform_points(ccrs.PlateCarree(),
                                  numpy.array(llcrnrlon),
                                  numpy.array(llcrnrlat))[0]
        ur = crs.transform_points(ccrs.PlateCarree(),
                                  numpy.array(urcrnrlon),
                                  numpy.array(urcrnrlat))[0]
    else:
        # one point that can be viewed from satellite
        lon0 = self.projection['satellite_lon'].get('degrees')
        lat0 = self.projection['satellite_lat'].get('degrees')
        # projected coordinates in both systems of this point
        pos_crs0 = crs.transform_points(ccrs.PlateCarree(),
                                        numpy.array(lon0),
                                        numpy.array(lat0))[0]
        pos_epy0 = self.ll2xy(lon0, lat0)
        # Corners
        ll = self.ij2xy(*self.gimme_corners_ij(subzone)['ll'])
        ur = self.ij2xy(*self.gimme_corners_ij(subzone)['ur'])
        ll = (ll[0] - pos_epy0[0] + pos_crs0[0],
              ll[1] - pos_epy0[1] + pos_crs0[1])
        ur = (ur[0] - pos_epy0[0] + pos_crs0[0],
              ur[1] - pos_epy0[1] + pos_crs0[1])
    return (ll[0], ur[0], ll[1], ur[1])