Python pyproj.CRS Examples

def _parse_crs(value):
        """ Method for parsing different inputs representing the same CRS enum. Examples:

        - 4326
        - 'EPSG:3857'
        - {'init': 32633}
        - pyproj.CRS(32743)
        """
        if isinstance(value, dict) and 'init' in value:
            value = value['init']
        if isinstance(value, pyproj.CRS):
            if value == CRSMeta._UNSUPPORTED_CRS:
                raise ValueError('sentinelhub-py supports only WGS 84 coordinate reference system with '
                                 'coordinate order lng-lat. However pyproj.CRS(4326) has coordinate order lat-lng')

            value = value.to_epsg()

        if isinstance(value, int):
            return str(value)
        if isinstance(value, str):
            return value.lower().strip('epsg: ')
        return value 

def test_proj4_radius_parameters_provided(self):
        """Test proj4_radius_parameters with a/b."""
        from pyresample import utils
        a, b = utils.proj4.proj4_radius_parameters(
            '+proj=stere +a=6378273 +b=6356889.44891',
        )
        np.testing.assert_almost_equal(a, 6378273)
        np.testing.assert_almost_equal(b, 6356889.44891)

        # test again but force pyproj <2 behavior
        with mock.patch.object(utils.proj4, 'CRS', None):
            a, b = utils.proj4.proj4_radius_parameters(
                '+proj=stere +a=6378273 +b=6356889.44891',
            )
            np.testing.assert_almost_equal(a, 6378273)
            np.testing.assert_almost_equal(b, 6356889.44891) 

def pyproj_transform(x, y, in_crs, out_crs, z=None):
    """
    Wrapper around pyproj to convert coordinates from an EPSG system to another.

    Args:
        x (scalar or array): x coordinate(s), expressed in in_crs
        y (scalar or array): y coordinate(s), expressed in in_crs
        in_crs (pyproj.crs.CRS or int): input coordinate reference system or EPSG code
        out_crs (pyproj.crs.CRS or int): output coordinate reference system or EPSG code
        z (scalar or array): z coordinate(s), expressed in in_crs

    Returns:
        scalar or array: x coordinate(s), expressed in out_crs
        scalar or array: y coordinate(s), expressed in out_crs
        scalar or array (optional if z): z coordinate(s), expressed in out_crs
    """
    transformer = pyproj.Transformer.from_crs(in_crs, out_crs, always_xy=True)
    if z is None:
        return transformer.transform(x, y)
    else:
        return transformer.transform(x, y, z) 

def pyproj_crs(projparams):
    """
    Wrapper around pyproj to return a pyproj.crs.CRS object that corresponds
    to the given parameters

    Args:
        projparams (int, str, dict): CRS parameters

    Returns:
        pyproj.crs.CRS: object that defines a CRS
    """
    if isinstance(projparams, str):
        try:
            projparams = int(projparams)
        except (ValueError, TypeError):
            pass
    return pyproj.crs.CRS(projparams) 

def rasterio_crs(projparams):
    """
    Return a rasterio.crs.CRS object that corresponds to the given parameters.
    See: https://pyproj4.github.io/pyproj/stable/crs_compatibility.html#converting-from-pyproj-crs-crs-to-rasterio-crs-crs

    Args:
        projparams (int, str, dict, pyproj.CRS): PROJ parameters

    Returns:
        rasterio.crs.CRS: object that can be used with rasterio
    """
    proj_crs = pyproj_crs(projparams)
    if LooseVersion(rasterio.__gdal_version__) < LooseVersion("3.0.0"):
        rio_crs = RioCRS.from_wkt(proj_crs.to_wkt(WktVersion.WKT1_GDAL))
    else:
        rio_crs = RioCRS.from_wkt(proj_crs.to_wkt())
    return rio_crs 

def test_crs_parsing(self):
        test_cases = [
            (4326, CRS.WGS84),
            ('4326', CRS.WGS84),
            ('EPSG:3857', CRS.POP_WEB),
            ({'init': 'EPSG:32638'}, CRS.UTM_38N),
            (pyproj.CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'), CRS.WGS84),
            (pyproj.CRS(3857), CRS.POP_WEB),
        ]
        for parse_value, expected_result in test_cases:
            with self.subTest(msg=str(parse_value)):
                parsed_result = CRS(parse_value)
                self.assertEqual(parsed_result, expected_result)

        with self.assertRaises(ValueError):
            CRS(pyproj.CRS(4326)) 

def test_utm(self):
        known_values = (
            (13, 46, '32633'),
            (13, 0, '32633'),
            (13, -45, '32733'),
            (13, 0, '32633'),
            (13, -0.0001, '32733'),
            (13, -46, '32733')
        )
        for known_val in known_values:
            lng, lat, epsg = known_val
            with self.subTest(msg=epsg):
                crs = CRS.get_utm_from_wgs84(lng, lat)
                self.assertEqual(epsg, crs.value,
                                 msg="Expected {}, got {} for lng={},lat={}".format(epsg, crs.value, str(lng),
                                                                                    str(lat))) 

def test_proj4_radius_parameters_spherical(self):
        """Test proj4_radius_parameters in case of a spherical earth."""
        from pyresample import utils
        a, b = utils.proj4.proj4_radius_parameters(
            '+proj=stere +R=6378273',
        )
        np.testing.assert_almost_equal(a, 6378273.)
        np.testing.assert_almost_equal(b, 6378273.)

        # test again but force pyproj <2 behavior
        with mock.patch.object(utils.proj4, 'CRS', None):
            a, b = utils.proj4.proj4_radius_parameters(
                '+proj=stere +R=6378273',
            )
            np.testing.assert_almost_equal(a, 6378273.)
            np.testing.assert_almost_equal(b, 6378273.) 

def test_proj4_radius_parameters_default(self):
        """Test proj4_radius_parameters with default parameters."""
        from pyresample import utils
        a, b = utils.proj4.proj4_radius_parameters(
            '+proj=lcc +lat_0=10 +lat_1=10',
        )
        # WGS84
        np.testing.assert_almost_equal(a, 6378137.)
        np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)

        # test again but force pyproj <2 behavior
        with mock.patch.object(utils.proj4, 'CRS', None):
            a, b = utils.proj4.proj4_radius_parameters(
                '+proj=lcc +lat_0=10 +lat_1=10',
            )
            # WGS84
            np.testing.assert_almost_equal(a, 6378137.)
            np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) 

def test_proj4_radius_parameters_ellps(self):
        """Test proj4_radius_parameters with ellps."""
        from pyresample import utils
        a, b = utils.proj4.proj4_radius_parameters(
            '+proj=stere +ellps=WGS84',
        )
        np.testing.assert_almost_equal(a, 6378137.)
        np.testing.assert_almost_equal(b, 6356752.314245, decimal=6)

        # test again but force pyproj <2 behavior
        with mock.patch.object(utils.proj4, 'CRS', None):
            a, b = utils.proj4.proj4_radius_parameters(
                '+proj=stere +ellps=WGS84',
            )
            np.testing.assert_almost_equal(a, 6378137.)
            np.testing.assert_almost_equal(b, 6356752.314245, decimal=6) 

def test_get_area_def_from_raster_extracts_proj_id(self):
        from rasterio.crs import CRS
        from pyresample import utils
        crs = CRS(init='epsg:3857')
        source = tmptiff(crs=crs)
        area_def = utils.rasterio.get_area_def_from_raster(source)
        self.assertEqual(area_def.proj_id, 'WGS 84 / Pseudo-Mercator') 

def get_proj(prj_code):
    """
      Helper method for handling projection codes that are unknown to pyproj

      Args:
          prj_code (str): an epsg proj code

      Returns:
          projection: a pyproj projection
    """
    if prj_code in CUSTOM_PRJ:
        proj = pyproj.CRS(CUSTOM_PRJ[prj_code])
    else:
        proj = pyproj.CRS(prj_code)
    return proj 

def friendly_crs_equal(expected, actual, keys=None, use_obj=True, use_wkt=True):
    """Test if two projection definitions are equal.

    The main purpose of this function is to help manage differences
    between pyproj versions. Depending on the version installed and used
    pyresample may provide a different `proj_dict` or other similar
    CRS definition.

    Args:
        expected (dict, str, pyproj.crs.CRS): Expected CRS definition as
            a PROJ dictionary or string or CRS object.
        actual (dict, str, pyproj.crs.CRS): Actual CRS definition
        keys (list): Specific PROJ parameters to look for. Only takes effect
            if `use_obj` is `False`.
        use_obj (bool): Use pyproj's CRS object to test equivalence. Default
            is True.
        use_wkt (bool): Increase likely hood of making CRS objects equal by
            converting WellKnownText before converting to the final CRS
            object. Requires `use_obj`. Defaults to True.

    """
    if CRS is not None and use_obj:
        if hasattr(expected, 'crs'):
            expected = expected.crs
        if hasattr(actual, 'crs'):
            actual = actual.crs
        expected_crs = CRS(expected)
        actual_crs = CRS(actual)
        if use_wkt:
            expected_crs = CRS(expected_crs.to_wkt())
            actual_crs = CRS(actual_crs.to_wkt())
        return expected_crs == actual_crs
    raise NotImplementedError("""TODO""") 

def __call__(cls, crs_value, *args, **kwargs):
        """ This is executed whenever CRS('something') is called
        """
        # pylint: disable=signature-differs
        crs_value = cls._parse_crs(crs_value)

        if isinstance(crs_value, str) and not cls.has_value(crs_value) and crs_value.isdigit() and len(crs_value) >= 4:
            crs_name = 'EPSG_{}'.format(crs_value)
            extend_enum(cls, crs_name, crs_value)

        return super().__call__(crs_value, *args, **kwargs) 

def test_get_area_def_from_raster_non_georef_respects_proj_dict(self):
        from pyresample import utils
        from affine import Affine
        transform = Affine(300.0379266750948, 0.0, 101985.0,
                           0.0, -300.041782729805, 2826915.0)
        source = tmptiff(transform=transform)
        proj_dict = {'init': 'epsg:3857'}
        area_def = utils.rasterio.get_area_def_from_raster(source, proj_dict=proj_dict)
        if utils.is_pyproj2():
            from pyproj import CRS
            proj_dict = CRS(3857).to_dict()
        self.assertDictEqual(area_def.proj_dict, proj_dict) 

def crs_bbx(ll_poly, crs=None):
    """
    Compute the UTM bounding box of a given (lon, lat) polygon.

    Args:
        ll_poly ()
        crs (pyproj.crs.CRS): pyproj CRS object. If not specified, the default CRS of the UTM
            zone for the given geography is used.

    Returns:
       4-tuple with easting min/max and northing min/max
    """
    if not crs:
        utm_zone = compute_utm_zone(*ll_poly.mean(axis=0))
        epsg = epsg_code_from_utm_zone(utm_zone)
        crs = pyproj_crs(epsg)

    # convert lon lat polygon to target CRS
    easting, northing = pyproj.transform(pyproj.Proj(init="epsg:4326"),
                                         crs, ll_poly[:, 0], ll_poly[:, 1])

    # return UTM bounding box
    east_min = min(easting)
    east_max = max(easting)
    nort_min = min(northing)
    nort_max = max(northing)
    return east_min, east_max, nort_min, nort_max 

def trajectories_from_qgis_point_layer(layer, time_field_name, trajectory_id_field, time_format):
    names = [field.name() for field in layer.fields()]
    data = []
    for feature in layer.getFeatures():
        my_dict = {}
        for i, a in enumerate(feature.attributes()):
            if names[i] == time_field_name:
                if type(a) == QtCore.QDateTime:
                    my_dict[names[i]] = a.toPyDateTime()
                else:
                    my_dict[names[i]] = datetime.strptime(a, time_format)
            else:
                my_dict[names[i]] = a
        x = feature.geometry().asPoint().x()
        y = feature.geometry().asPoint().y()
        my_dict['geometry'] = Point((x, y))
        data.append(my_dict)
    df = pd.DataFrame(data).set_index(time_field_name)
    crs = CRS(int(layer.sourceCrs().geographicCrsAuthId().split(':')[1]))
    geo_df = GeoDataFrame(df, crs=crs)
    df_by_id = dict(tuple(geo_df.groupby(trajectory_id_field)))
    trajectories = []
    for key, value in df_by_id.items():
        traj = Trajectory(key, value)
        trajectories.append(traj)
    return trajectories 

def _check_crs(input_crs, return_rasterio=False):
    """Convert CRS to the ``pyproj.CRS`` object passed by ``solaris``."""
    if not isinstance(input_crs, pyproj.CRS) and input_crs is not None:
        out_crs = pyproj.CRS(input_crs)
    else:
        out_crs = input_crs

    if return_rasterio:
        if LooseVersion(rasterio.__gdal_version__) >= LooseVersion("3.0.0"):
            out_crs = rasterio.crs.CRS.from_wkt(out_crs.to_wkt())
        else:
            out_crs = rasterio.crs.CRS.from_wkt(out_crs.to_wkt("WKT1_GDAL"))

    return out_crs 

def test_custom_crs(self):
        for incorrect_value in ['string', -1, 999, None]:
            with self.assertRaises(ValueError):
                CRS(incorrect_value)

        for correct_value in [3035, 'EPSG:3035', 10000]:
            CRS(CRS(correct_value))

            new_enum_value = str(correct_value).lower().strip('epsg: ')
            self.assertTrue(CRS.has_value(new_enum_value)) 

def __new__(mcs, cls, bases, classdict):
        """ This is executed at the beginning of runtime when CRS class is created
        """
        for direction, direction_value in [('N', '6'), ('S', '7')]:
            for zone in range(1, 61):
                classdict['UTM_{}{}'.format(zone, direction)] = '32{}{}'.format(direction_value, str(zone).zfill(2))

        return super().__new__(mcs, cls, bases, classdict) 

def __str__(self):
        """ Method for casting CRS enum into string
        """
        return self.ogc_string() 

def __repr__(self):
        """ Method for retrieving CRS enum representation
        """
        return "CRS('{}')".format(self.value) 

def epsg(self):
        """ EPSG code property

        :return: EPSG code of given CRS
        :rtype: int
        """
        return int(self.value) 

def ogc_string(self):
        """ Returns a string of the form authority:id representing the CRS.

        :param self: An enum constant representing a coordinate reference system.
        :type self: CRS
        :return: A string representation of the CRS.
        :rtype: str
        """
        return 'EPSG:{}'.format(CRS(self).value) 

def opengis_string(self):
        """ Returns an URL to OGC webpage where the CRS is defined

        :return: An URL with CRS definition
        :rtype: str
        """
        return 'http://www.opengis.net/def/crs/EPSG/0/{}'.format(self.epsg) 

def is_utm(self):
        """ Checks if crs is one of the 64 possible UTM coordinate reference systems.

        :param self: An enum constant representing a coordinate reference system.
        :type self: CRS
        :return: `True` if crs is UTM and `False` otherwise
        :rtype: bool
        """
        return self.name.startswith('UTM') 

def projection(self):
        """ Returns a projection in form of pyproj class. For better time performance it will cache results of
        5 most recently used CRS classes.

        :return: pyproj projection class
        :rtype: pyproj.Proj
        """
        return pyproj.Proj(self._get_pyproj_projection_def(), preserve_units=True) 

def get_transform_function(self, other):
        """ Returns a function for transforming geometrical objects from one CRS to another. The function will support
        transformations between any objects that pyproj supports.
        For better time performance this method will cache results of 10 most recently used pairs of CRS classes.

        :param self: Initial CRS
        :type self: CRS
        :param other: Target CRS
        :type other: CRS
        :return: A projection function obtained from pyproj package
        :rtype: function
        """
        return pyproj.Transformer.from_proj(self.projection(), other.projection(), skip_equivalent=True).transform 

def get_utm_from_wgs84(lng, lat):
        """ Convert from WGS84 to UTM coordinate system

        :param lng: Longitude
        :type lng: float
        :param lat: Latitude
        :type lat: float
        :return: UTM coordinates
        :rtype: tuple
        """
        _, _, zone, _ = utm.from_latlon(lat, lng)
        direction = 'N' if lat >= 0 else 'S'
        return CRS['UTM_{}{}'.format(str(zone), direction)] 

def _get_pyproj_projection_def(self):
        """ Returns a pyproj crs definition

        For WGS 84 it ensures lng-lat order
        """
        return '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs' if self is CRS.WGS84 else self.ogc_string()