Repository: euro-cordex/py-cordex

Issue:
use more general grid mapping
We would like to have more general grid mappings. In this example, we use cartopy to compute the transformed coordinates in the global lat lon grid:

```python
import xarray as xr
from cartopy import crs as ccrs
import cordex as cx

eur11 = cx.cordex_domain('EUR-11', dummy='topo')
eur11

# reconstruct lat lon coordinates using cartopy
pole = ccrs.RotatedPole(-162.0, 39.25)
latlon = ccrs.PlateCarree()

lat_stack, lon_stack = xr.broadcast(eur11.rlat, eur11.rlon)
result = latlon.transform_points(pole, lon_stack.values, lat_stack.values)
lon, lat = result[:,:,0], result[:,:,1]
```
check differences to our function:
```python
np.abs(eur11.lon - lon).max(), np.abs(eur11.lat - lat).max()
```


Discussion:
codecov-commenter: # [Codecov](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex) Report
> Merging [#71](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex) (0af5de4) into [master](https://codecov.io/gh/euro-cordex/py-cordex/commit/11a037844bd19677ad6f2a5e63eb526db9f19988?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex) (11a0378) will **decrease** coverage by `4.36%`.
> The diff coverage is `74.54%`.

:mega: This organization is not using Codecov’s [GitHub App Integration](https://github.com/apps/codecov). We recommend you install it so Codecov can continue to function properly for your repositories. [Learn more](https://about.codecov.io/blog/codecov-is-updating-its-github-integration/?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex)

<details><summary>Additional details and impacted files</summary>


[![Impacted file tree graph](https://codecov.io/gh/euro-cordex/py-cordex/pull/71/graphs/tree.svg?width=650&height=150&src=pr&token=D6BAMOZZT3&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex)](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex)

```diff
@@            Coverage Diff             @@
##           master      #71      +/-   ##
==========================================
- Coverage   77.85%   73.49%   -4.36%     
==========================================
  Files          20       21       +1     
  Lines        1016     1030      +14     
==========================================
- Hits          791      757      -34     
- Misses        225      273      +48     
```

| Flag | Coverage Δ | |
|---|---|---|
| unittests | `73.49% <74.54%> (-4.36%)` | :arrow_down: |

Flags with carried forward coverage won't be shown. [Click here](https://docs.codecov.io/docs/carryforward-flags?utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex#carryforward-flags-in-the-pull-request-comment) to find out more.

| [Impacted Files](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex) | Coverage Δ | |
|---|---|---|
| [cordex/core/utils.py](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex#diff-Y29yZGV4L2NvcmUvdXRpbHMucHk=) | `50.00% <ø> (-27.42%)` | :arrow_down: |
| [cordex/core/domain.py](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex#diff-Y29yZGV4L2NvcmUvZG9tYWluLnB5) | `69.36% <36.36%> (-5.86%)` | :arrow_down: |
| [cordex/\_\_init\_\_.py](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex#diff-Y29yZGV4L19faW5pdF9fLnB5) | `80.00% <100.00%> (ø)` | |
| [cordex/core/transform.py](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex#diff-Y29yZGV4L2NvcmUvdHJhbnNmb3JtLnB5) | `100.00% <100.00%> (ø)` | |
| [cordex/preprocessing/preprocessing.py](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex#diff-Y29yZGV4L3ByZXByb2Nlc3NpbmcvcHJlcHJvY2Vzc2luZy5weQ==) | `63.03% <0.00%> (-17.07%)` | :arrow_down: |
| [cordex/cmor/cmor.py](https://codecov.io/gh/euro-cordex/py-cordex/pull/71?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=euro-cordex#diff-Y29yZGV4L2Ntb3IvY21vci5weQ==) | `72.35% <0.00%> (+0.33%)` | :arrow_up: |

</details>

Patch:
diff --git a/.github/workflows/ci.yaml b/.github/workflows/ci.yaml
index 33b9d58..ff56619 100644
--- a/.github/workflows/ci.yaml
+++ b/.github/workflows/ci.yaml
@@ -26,7 +26,7 @@ jobs:
       matrix:
         os: ["ubuntu-latest"]
         # python versions
-        python-version: ["3.7", "3.8", "3.9", "3.10", "3.11"]
+        python-version: ["3.8", "3.9", "3.10", "3.11"]
     steps:
       - uses: actions/checkout@v3
         with:
diff --git a/ci/requirements/docs.yml b/ci/requirements/docs.yml
index ea29166..0f1e871 100644
--- a/ci/requirements/docs.yml
+++ b/ci/requirements/docs.yml
@@ -18,6 +18,7 @@ dependencies:
   - ipython
   - cftime
   - pooch
+  - pyproj
   - sphinx==4.2.0
   - sphinx_rtd_theme==1.0.0
   - sphinx-book-theme
diff --git a/ci/requirements/environment.yml b/ci/requirements/environment.yml
index 36fc0d3..f509d4e 100644
--- a/ci/requirements/environment.yml
+++ b/ci/requirements/environment.yml
@@ -4,12 +4,14 @@ channels:
   - nodefaults
 dependencies:
   - cartopy
+  - pyproj>=3.3.0 # rotated pole bug: https://github.com/pyproj4/pyproj/issues/948
   - matplotlib-base
   - numpy<1.24 # fixed until next cmor release
   - pooch
-  - regionmask
+  - regionmask>=0.9.0
   - setuptools
   - xarray
+  - cf_xarray>=0.8.0
   - cftime
   - netCDF4
   - xesmf
diff --git a/cordex/__init__.py b/cordex/__init__.py
index 79eb0a7..83d529e 100644
--- a/cordex/__init__.py
+++ b/cordex/__init__.py
@@ -1,22 +1,15 @@
-# -*- coding: utf-8 -*-
-
 import pkg_resources
 
-from . import core, regions, tables, tutorial  # , cmor, preprocessing
-from .core.domain import (
-    cordex_domain,
-    create_dataset,
-    domain_info,
+from . import core, regions, tables, tutorial
+from .core.domain import cordex_domain, create_dataset, domain_info, vertices
+from .core.transform import (
+    map_crs,
     rotated_coord_transform,
-    vertices,
+    transform,
+    transform_coords,
 )
-from .core.utils import map_crs
 from .tables import domains, ecmwf
 
-# from .version import version
-
-# __version__ = version
-
 try:
     __version__ = pkg_resources.get_distribution("py-cordex").version
 except Exception:
@@ -36,6 +29,8 @@ __all__ = [
     "rotated_coord_transform",
     "vertices",
     "map_crs",
+    "transform",
+    "transform_coords",
     "domains",
     "ecmwf",
 ]
diff --git a/cordex/core/domain.py b/cordex/core/domain.py
index 5ab4766..de8fa31 100644
--- a/cordex/core/domain.py
+++ b/cordex/core/domain.py
@@ -19,13 +19,16 @@ Example:
 
 from warnings import warn
 
+import cf_xarray as cfxr
 import numpy as np
 import pandas as pd
 import xarray as xr
+from pyproj import CRS
 
 from ..tables import domains
 from . import cf, utils
 from .config import nround
+from .transform import grid_mapping, transform, transform_coords
 
 
 def domain_names(table_name=None):
@@ -62,7 +65,7 @@ def cordex_domain(
     Parameters
     ----------
     short_name:
-        Name of the Cordex Domain.
+        Domain identifier.
     dummy : str or logical
         Name of dummy field, if dummy=topo, the cdo topo operator will be
         used to create some dummy topography data. dummy data is useful for
@@ -181,18 +184,16 @@ def create_dataset(
     except:
         pass
 
-    x, y = _init_grid(nlon, nlat, dlon, dlat, ll_lon, ll_lat)
+    x, y = _lin_coord(nlon, dlon, ll_lon), _lin_coord(nlat, dlat, ll_lat)
+
     if rotated is True:
-        lon, lat = rotated_coord_transform(*_stack(x, y), pollon, pollat)
-        pole = _grid_mapping(pollon, pollat)
-        if mapping_name is None:
-            mapping_name = cf.DEFAULT_MAPPING_NCVAR
         return _get_rotated_dataset(
             x,
             y,
-            lon,
-            lat,
-            pole,
+            #  lon,
+            #  lat,
+            pollon,
+            pollat,
             add_vertices=add_vertices,
             dummy=dummy,
             mapping_name=mapping_name,
@@ -266,124 +267,89 @@ def _get_regular_dataset(
             ds[mapping_name].grid_north_pole_longitude,
             ds[mapping_name].grid_north_pole_latitude,
         )
+        # v = vertices(ds, ccrs.RotatedPole(*pole))
         v = vertices(ds.rlon, ds.rlat, ccrs.RotatedPole(*pole))
         ds = xr.merge([ds, v])
         ds[cf.LON_NAME].attrs["bounds"] = cf.LON_BOUNDS
         ds[cf.LAT_NAME].attrs["bounds"] = cf.LAT_BOUNDS
 
-    if dummy:
-        if dummy is True:
-            dummy_name = "dummy"
-        else:
-            dummy_name = dummy
-        dummy = xr.DataArray(
-            data=np.zeros((ds.lat.size, ds.lon.size)),
-            coords=(ds.lat, ds.lon),
-        )
-        dummy.attrs = {"coordinates": "lat lon"}
-        ds[dummy_name] = dummy
-        if dummy_name == "topo":
-            # use cdo to create dummy topography data.
-            from cdo import Cdo
-
-            tmp = utils.get_tempfile()
-            ds.to_netcdf(tmp)
-            topo = Cdo().topo(tmp, returnXDataset=True)["topo"]
-            ds[dummy_name] = xr.DataArray(
-                data=topo.values,
-                coords=(ds.lat, ds.lon),
-            )
-            ds[dummy_name].attrs.update(topo.attrs)
+    if dummy is not None:
+        ds = _add_dummy(ds, dummy)
     return ds
 
 
 def _get_rotated_dataset(
     rlon,
     rlat,
-    lon,
-    lat,
-    pole,
+    pollon,
+    pollat,
     add_vertices=False,
     dummy=None,
     mapping_name=None,
     attrs=None,
 ):
-    data_vars = {mapping_name: pole}
+    mapping = grid_mapping(pollon, pollat, mapping_name)
+    # lon, lat = transform(x, y, src_crs=CRS.from_cf(mapping.attrs))
 
     ds = xr.Dataset(
-        data_vars=data_vars,
+        data_vars={mapping.name: mapping},
         coords=dict(
             rlon=(cf.RLON_NAME, rlon),
             rlat=([cf.RLAT_NAME], rlat),
-            lon=([cf.RLAT_NAME, cf.RLON_NAME], lon),
-            lat=([cf.RLAT_NAME, cf.RLON_NAME], lat),
         ),
         attrs=attrs,
     )
 
-    for key, coord in ds.coords.items():
-        coord.encoding["_FillValue"] = None
-        coord.attrs = cf.coords[key]
+    lon, lat = transform(ds.rlon, ds.rlat, src_crs=CRS.from_cf(mapping.attrs))
+    ds = ds.assign_coords(lon=lon, lat=lat)
 
     if add_vertices is True:
-        from cartopy import crs as ccrs
-
-        pole = (
-            ds[mapping_name].grid_north_pole_longitude,
-            ds[mapping_name].grid_north_pole_latitude,
-        )
-        v = vertices(ds.rlon, ds.rlat, ccrs.RotatedPole(*pole))
+        v = vertices(ds.rlon, ds.rlat, CRS.from_cf(mapping.attrs))
         ds = xr.merge([ds, v])
         ds[cf.LON_NAME].attrs["bounds"] = cf.LON_BOUNDS
         ds[cf.LAT_NAME].attrs["bounds"] = cf.LAT_BOUNDS
 
-    if dummy:
-        if dummy is True:
-            dummy_name = "dummy"
-        else:
-            dummy_name = dummy
-        dummy = xr.DataArray(
-            data=np.zeros((ds.rlat.size, ds.rlon.size)),
-            coords=(ds.rlat, ds.rlon),
-        )
-        dummy.attrs = {"grid_mapping": mapping_name, "coordinates": "lat lon"}
-        ds[dummy_name] = dummy
-        if dummy_name == "topo":
-            # use cdo to create dummy topography data.
-            from cdo import Cdo
-
-            tmp = utils.get_tempfile()
-            ds.to_netcdf(tmp)
-            topo = Cdo().topo(tmp, returnXDataset=True)["topo"]
-            ds[dummy_name] = xr.DataArray(
-                data=topo.values,
-                coords=(ds.rlat, ds.rlon),
-            )
-            ds[dummy_name].attrs.update(topo.attrs)
+    for key, coord in ds.coords.items():
+        coord.encoding["_FillValue"] = None
+        coord.attrs = cf.coords[key]
+
+    if dummy is not None:
+        ds = _add_dummy(ds, dummy)
     return ds
 
 
-def _grid_mapping(pollon, pollat, mapping_name=None):
-    """creates a grid mapping DataArray object"""
-    if mapping_name is None:
-        mapping_name = cf.DEFAULT_MAPPING_NCVAR
-    da = xr.DataArray(np.zeros((), dtype=np.int32))
-    attrs = cf.mapping.copy()
-    attrs["grid_north_pole_longitude"] = pollon
-    attrs["grid_north_pole_latitude"] = pollat
-    da.attrs = attrs
-    return da
+def _add_dummy(ds, name=True):
+    if name is True:
+        name = "dummy"
+
+    attrs = {"coordinates": "lat lon"}
+    try:
+        attrs["grid_mapping"] = ds.cf["grid_mapping"].name
+    except KeyError:
+        pass
+    # this is required for CDO to understand coordinates
+    ds[name] = xr.DataArray(
+        data=np.zeros((ds.cf.dims["Y"], ds.cf.dims["X"])),
+        coords=(ds.cf["Y"], ds.cf["X"]),
+        attrs=attrs,
+    )
+    if name == "topo":
+        # use cdo to create dummy topography data.
+        from cdo import Cdo
+
+        tmp = utils.get_tempfile()
+        ds.to_netcdf(tmp)
+        topo = Cdo().topo(tmp, returnXDataset=True)["topo"]
+        ds[name] = xr.DataArray(
+            data=topo.values, coords=(ds.cf["Y"], ds.cf["X"]), attrs=attrs
+        )
+        ds[name].attrs.update(topo.attrs)
+    return ds
 
 
-def _init_grid(nlon, nlat, dlon, dlat, ll_lon, ll_lat):
+def _lin_coord(nx, dx, x0, dtype=np.float64):
     """create coordinate arrays from lower left longitude and latitude"""
-    rlon = np.array(
-        [round(ll_lon + i * dlon, nround) for i in range(0, nlon)], dtype=np.float64
-    )
-    rlat = np.array(
-        [round(ll_lat + i * dlat, nround) for i in range(0, nlat)], dtype=np.float64
-    )
-    return rlon, rlat
+    return np.round(np.arange(0, nx, dtype=dtype) * dx + x0, nround)
 
 
 def _stack(x, y):
@@ -395,178 +361,6 @@ def _stack(x, y):
     )
 
 
-def rotated_coord_transform(lon, lat, np_lon, np_lat, direction="rot2geo"):
-    """Transforms a coordinate into a rotated grid coordinate and vice versa.
-
-    The coordinates have to be given in degree and will be returned in degree.
-
-    Parameters
-    ----------
-    lon : float array like
-        Longitude coordinate.
-    lat : float array like
-        Latitude coordinate.
-    np_lon : float array like
-        Longitude coordinate of the rotated north pole.
-    np_lat : float array like
-        Latitude coordinate of the rotated north pole.
-    direction : str
-        Direction of the rotation.
-        Options are: 'rot2geo' (default) for a transformation to regular
-        coordinates from rotated. 'geo2rot' transforms regular coordinates
-        to rotated.
-
-    Returns
-    -------
-    lon_new : array like
-        New longitude coordinate.
-    lat_new : array like
-        New latitude coordinate.
-    """
-    # Convert degrees to radians
-    lon = np.deg2rad(lon)
-    lat = np.deg2rad(lat)
-
-    theta = 90.0 - np_lat  # Rotation around y-axis
-    phi = np_lon + 180.0  # Rotation around z-axis
-
-    # Convert degrees to radians
-    phi = np.deg2rad(phi)
-    theta = np.deg2rad(theta)
-
-    # Convert from spherical to cartesian coordinates
-    x = np.cos(lon) * np.cos(lat)
-    y = np.sin(lon) * np.cos(lat)
-    z = np.sin(lat)
-
-    # Regular -> Rotated
-    if direction == "geo2rot":
-        x_new = (
-            np.cos(theta) * np.cos(phi) * x
-            + np.cos(theta) * np.sin(phi) * y
-            + np.sin(theta) * z
-        )
-        y_new = -np.sin(phi) * x + np.cos(phi) * y
-        z_new = (
-            -np.sin(theta) * np.cos(phi) * x
-            - np.sin(theta) * np.sin(phi) * y
-            + np.cos(theta) * z
-        )
-
-    # Rotated -> Regular
-    elif direction == "rot2geo":
-        phi = -phi
-        theta = -theta
-
-        x_new = (
-            np.cos(theta) * np.cos(phi) * x
-            + np.sin(phi) * y
-            + np.sin(theta) * np.cos(phi) * z
-        )
-        y_new = (
-            -np.cos(theta) * np.sin(phi) * x
-            + np.cos(phi) * y
-            - np.sin(theta) * np.sin(phi) * z
-        )
-        z_new = -np.sin(theta) * x + np.cos(theta) * z
-
-    # Convert cartesian back to spherical coordinates
-    lon_new = np.arctan2(y_new, x_new)
-    lat_new = np.arcsin(z_new)
-
-    # Convert radians back to degrees
-    lon_new = np.rad2deg(lon_new)
-    lat_new = np.rad2deg(lat_new)
-
-    return lon_new, lat_new
-
-
-def _map_crs(lon, lat, src_crs, trg_crs=None):
-    """coordinate transformation of longitude and latitude
-
-    Transforms the coordinates lat, lon from the transform crs
-    into the projection crs using cartopy.crs.
-
-    Parameters
-    ----------
-    lon : float array like
-        Longitude coordinate.
-    lat : float array like
-        Latitude coordinate.
-    src_crs : cartopy.crs
-        Source coordinate reference system into which lat and lon
-        should are defined.
-    trg_crs : cartopy.crs
-        Target coordinate reference system in which lat and lon
-        should be transformed. If `None`, `PlateCarree` is used.
-
-    Returns
-    -------
-    lon : array like
-        Projected longitude coordinate.
-    lat : array like
-        Projected latitude coordinate.
-
-    """
-
-    from cartopy import crs as ccrs
-
-    if trg_crs is None:
-        trg_crs = ccrs.PlateCarree()
-    # latlon = ccrs.PlateCarree()
-    # lon_stack, lat_stack = xr.broadcast(lon, lat)
-    lon_stack = np.broadcast_to(lon, (lat.shape[0], lon.shape[0])).T
-    lat_stack = np.broadcast_to(lat, (lon.shape[0], lat.shape[0]))
-    result = trg_crs.transform_points(src_crs, lon_stack, lat_stack)
-    return result[:, :, 0], result[:, :, 1]
-
-
-# wrapper function for xarray.apply_ufunc
-def map_crs(lon, lat, src_crs, trg_crs=None):
-    """coordinate transformation of longitude and latitude
-
-    Transforms the coordinates lat, lon from the transform crs
-    into the projection crs using cartopy.crs.
-
-    Parameters
-    ----------
-    lon : float array like
-        Longitude coordinate.
-    lat : float array like
-        Latitude coordinate.
-    src_crs : cartopy.crs
-        Source coordinate reference system into which lat and lon
-        should are defined.
-    trg_crs : cartopy.crs
-        Target coordinate reference system in which lat and lon
-        should be transformed. If `None`, `PlateCarree` is used.
-
-    Returns
-    -------
-    lon : xr.DataArray
-        Projected longitude coordinate with dims (lat, lon).
-    lat : xr.DataArray
-        Projected latitude coordinate with dims (lat, lon).
-
-    """
-    input_core_dims = [[lon.dims[0]], [lat.dims[0]]] + [[], []]
-    output_core_dims = 2 * [[lon.dims[0], lat.dims[0]]]
-    result = xr.apply_ufunc(
-        _map_crs,  # first the function
-        lon,  # now arguments in the order expected by 'interp1_np'
-        lat,
-        src_crs,
-        trg_crs,
-        input_core_dims=input_core_dims,  # list with one entry per arg
-        # [["rlat", "rlon"], ["rlat", "rlon"]],
-        output_core_dims=output_core_dims
-        # exclude_dims=set(("lat",)),  # dimensions allowed to change size. Must be set!
-    )
-    result[0].name = cf.LON_NAME
-    result[1].name = cf.LAT_NAME
-    return result
-
-
 def _dcoord(coord, include="left"):
     dcoord = coord.values[1:] - coord.values[:-1]
     if include in ["left", "both"]:
@@ -606,7 +400,7 @@ def bounds(coords):
     # return xr.merge([_bounds(coord).left.rename({"left": coord.name+"bounds"}) for coord in coords])
 
 
-def vertices(rlon, rlat, src_crs, trg_crs=None):
+def vertices_new(ds, src_crs, trg_crs=None):
     """Compute lon and lat vertices.
 
     Transformation of rlon vertices and rlat vertices
@@ -629,14 +423,27 @@ def vertices(rlon, rlat, src_crs, trg_crs=None):
     warn(
         "Order of vertices has changed since v0.3.2 to CF Conventions, see https://github.com/euro-cordex/py-cordex/issues/34"
     )
-    rlon_bounds = _bounds(rlon)
-    rlat_bounds = _bounds(rlat)
+    # ds = xr.merge([rlon_bounds, rlat_bounds])
+    rlon_bounds = ds.cf.add_bounds("rlon").rlon_bounds.drop(
+        "rlon_bounds"
+    )  # _bounds(ds.rlon)
+    rlat_bounds = ds.cf.add_bounds("rlat").rlat_bounds.drop(
+        "rlat_bounds"
+    )  # _bounds(ds.rlat)
     # maps each vertex to lat lon coordinates
     # order is counterclockwise starting from lower left vertex
-    v1 = map_crs(rlon_bounds.left, rlat_bounds.left, src_crs, trg_crs)
-    v2 = map_crs(rlon_bounds.right, rlat_bounds.left, src_crs, trg_crs)
-    v3 = map_crs(rlon_bounds.right, rlat_bounds.right, src_crs, trg_crs)
-    v4 = map_crs(rlon_bounds.left, rlat_bounds.right, src_crs, trg_crs)
+    v1 = transform(
+        rlon_bounds.isel(bounds=0), rlat_bounds.isel(bounds=0), src_crs, trg_crs
+    )
+    v2 = transform(
+        rlon_bounds.isel(bounds=1), rlat_bounds.isel(bounds=0), src_crs, trg_crs
+    )
+    v3 = transform(
+        rlon_bounds.isel(bounds=1), rlat_bounds.isel(bounds=1), src_crs, trg_crs
+    )
+    v4 = transform(
+        rlon_bounds.isel(bounds=0), rlat_bounds.isel(bounds=1), src_crs, trg_crs
+    )
     lon_vertices = xr.concat(
         [v1[0], v2[0], v3[0], v4[0]], dim=cf.BOUNDS_DIM
     ).transpose()
@@ -654,8 +461,8 @@ def vertices(rlon, rlat, src_crs, trg_crs=None):
     return xr.merge([lat_vertices, lon_vertices])
 
 
-def lon_lat_bounds_coordinates(ds, src_crs=None, trg_crs=None):
-    """Compute lon and lat bounds coordinates.
+def vertices(rlon, rlat, src_crs, trg_crs=None):
+    """Compute lon and lat vertices.
 
     Transformation of rlon vertices and rlat vertices
     into the target crs according to
@@ -674,32 +481,26 @@ def lon_lat_bounds_coordinates(ds, src_crs=None, trg_crs=None):
         lon_vertices and lat_vertices in target crs.
 
     """
-    if src_crs is None:
-        src_crs = get_pole_crs(ds)
-    rot_bounds = bounds_coordinates(ds, ("rlon", "rlat"))
-    lon_bounds, lat_bounds = map_crs(
-        rot_bounds.coords["rlon_b"], rot_bounds.coords["rlat_b"], src_crs, trg_crs
-    )
-    return ds.assign_coords(lon_b=lon_bounds.transpose(), lat_b=lat_bounds.transpose())
-
-
-def get_pole_crs(ds, grid_mapping_name="rotated_latitude_longitude"):
-    """Create cartopy crs from rotated pole dataset."""
-    from cartopy import crs as ccrs
-
-    pole = (
-        ds[grid_mapping_name].grid_north_pole_longitude,
-        ds[grid_mapping_name].grid_north_pole_latitude,
-    )
-    return ccrs.RotatedPole(*pole)
-
-
-def extend_coordinate(coord, n=1):
-    data = coord
-    # insert = np.array([data[0] - n*(data[1]-data[0]) + i*(data[1]-data[0]) for i in range(0, n)])
-    insert = np.array([data[0] - (n - i) * (data[1] - data[0]) for i in range(0, n)])
-    return insert
-    data_extended = np.insert(data, 0, data[0] - (data[1] - data[0]))
-    # return data_extended
-    # data_extended = np.append(data_extended, data[-1]+(data[-1]-data[-2]))
-    return xr.DataArray(data_extended, dims=coord.dims, name=coord.name)
+    rlon_bounds = _bounds(rlon)
+    rlat_bounds = _bounds(rlat)
+    # maps each vertex to lat lon coordinates
+    # order is counterclockwise starting from lower left vertex
+    v1 = transform(rlon_bounds.left, rlat_bounds.left, src_crs, trg_crs)
+    v2 = transform(rlon_bounds.right, rlat_bounds.left, src_crs, trg_crs)
+    v3 = transform(rlon_bounds.right, rlat_bounds.right, src_crs, trg_crs)
+    v4 = transform(rlon_bounds.left, rlat_bounds.right, src_crs, trg_crs)
+    lon_vertices = xr.concat(
+        [v1[0].T, v2[0].T, v3[0].T, v4[0].T], dim=cf.BOUNDS_DIM
+    ).transpose()
+    #    ..., "vertices"
+    # )
+    lat_vertices = xr.concat(
+        [v1[1].T, v2[1].T, v3[1].T, v4[1].T], dim=cf.BOUNDS_DIM
+    ).transpose()
+    #    ..., "vertices"
+    # )
+    lon_vertices.name = cf.LON_BOUNDS
+    lat_vertices.name = cf.LAT_BOUNDS
+    lon_vertices.attrs = cf.coords[cf.LON_BOUNDS]
+    lat_vertices.attrs = cf.coords[cf.LAT_BOUNDS]
+    return xr.merge([lat_vertices, lon_vertices])
diff --git a/cordex/core/transform.py b/cordex/core/transform.py
new file mode 100644
index 0000000..04f1301
--- /dev/null
+++ b/cordex/core/transform.py
@@ -0,0 +1,278 @@
+from warnings import warn
+
+import numpy as np
+import xarray as xr
+from pyproj import CRS, Transformer
+
+from . import cf
+
+xr.set_options(keep_attrs=True)
+
+
+def _map_crs(x_stack, y_stack, src_crs, trg_crs=None):
+    """coordinate transformation of longitude and latitude"""
+
+    from cartopy import crs as ccrs
+
+    if trg_crs is None:
+        trg_crs = ccrs.PlateCarree()
+    result = trg_crs.transform_points(src_crs, x_stack, y_stack)
+    return result[:, :, 0], result[:, :, 1]
+
+
+# wrapper function for xarray.apply_ufunc
+def map_crs(x, y, src_crs, trg_crs=None):
+    """coordinate transformation using cartopy
+
+    Transforms the coordinates x, y from the source crs
+    into the target crs using cartopy.
+
+    Parameters
+    ----------
+    x : float array like
+        x coordinate of source crs.
+    y : float array like
+        y coordinate of source crs.
+    src_crs : cartopy.crs
+        Source coordinate reference system in which x and y
+        are defined.
+    trg_crs : cartopy.crs
+        Target coordinate reference system into which x and y
+        should be transformed. If `None`, `PlateCarree` is used.
+
+    Returns
+    -------
+    x_map : xr.DataArray
+        Projected x coordinate.
+    y_map : xr.DataArray
+        Projected y coordinate.
+
+    """
+    warn(
+        "map_crs is deprecated, please use transform_xy instead",
+        DeprecationWarning,
+        stacklevel=2,
+    )
+    y_stack, x_stack = xr.broadcast(y, x)
+    input_core_dims = 2 * [list(x_stack.dims)] + [[], []]
+    output_core_dims = 2 * [list(x_stack.dims)]
+
+    result = xr.apply_ufunc(
+        _map_crs,  # first the function
+        x_stack,  # now arguments in the order expected by 'interp1_np'
+        y_stack,
+        src_crs,
+        trg_crs,
+        input_core_dims=input_core_dims,  # list with one entry per arg
+        # [["rlat", "rlon"], ["rlat", "rlon"]],
+        output_core_dims=output_core_dims
+        # exclude_dims=set(("lat",)),  # dimensions allowed to change size. Must be set!
+    )
+
+    result[0].name = "x_map"
+    result[1].name = "y_map"
+
+    return result
+
+
+def _transform(x, y, src_crs, trg_crs):
+    """helper function for transforming coordinates"""
+    # always_xy=True
+    # https://proj.org/faq.html#why-is-the-axis-ordering-in-proj-not-consistent
+    transformer = Transformer.from_crs(src_crs, trg_crs, always_xy=True)
+    xt, yt = transformer.transform(x, y)
+    return xt, yt
+
+
+def transform(x, y, src_crs, trg_crs=None):
+    """Coordinate transformation using pyproj.
+
+    Transforms the coordinates x, y from the source crs
+    into a target crs using pyproj.
+
+    Parameters
+    ----------
+    x : DataArray
+        X coordinate.
+    y : DataArray
+        Y coordinate.
+    src_crs : pyproj.CRS
+        Source coordinate reference system in which x and y are defined.
+    trg_crs : pyproj.CRS
+        Target coordinate reference system into which x and y
+        should be transformed. If not supplied, ``EPSG:4326`` is the default.
+
+    Returns
+    -------
+    xt : DataArray
+        Transformed x coordinate.
+    yt : DataArray
+        Transformed y coordinate.
+
+    """
+    if trg_crs is None:
+        # default target crs
+        trg_crs = CRS("EPSG:4326")
+    y_stack, x_stack = xr.broadcast(y, x)
+    input_core_dims = [x_stack.dims, y_stack.dims] + [[], []]
+    output_core_dims = [x_stack.dims, y_stack.dims]
+
+    xt, yt = xr.apply_ufunc(
+        _transform,
+        x_stack,
+        y_stack,
+        src_crs,
+        trg_crs,
+        input_core_dims=input_core_dims,
+        output_core_dims=output_core_dims,
+    )
+
+    xt.name = "xt"
+    yt.name = "yt"
+    xt.attrs = {"epsg": trg_crs.to_epsg()}
+    yt.attrs = {"epsg": trg_crs.to_epsg()}
+
+    return xt, yt
+
+
+def transform_coords(ds, src_crs=None, trg_crs=None, trg_dims=None):
+    """Transform X and Y coordinates of a Dataset.
+
+    The transformed coordinates will be added to the Dataset.
+
+    Parameters
+    ----------
+    ds : Dataset or DataArray
+        Dataset with input grid.
+    src_crs : pyproj.CRS
+        Source coordinate reference system in which X and Y are defined.
+        If not supplied, a `grid_mapping` variable should be available
+        to define the source CRS.
+    trg_crs : pyproj.CRS
+        Target coordinate reference system into which x and y
+        should be transformed. If not supplied, ``EPSG:4326`` is the default.
+    trg_dims: list or set
+        Names of the output coordinates.
+
+    Returns
+    -------
+    ds : Dataset or DataArray
+        Dataset with transformed coordinates.
+
+    """
+
+    if trg_crs is None:
+        # default target crs
+        trg_crs = CRS("EPSG:4326")
+    if trg_dims is None:
+        trg_dims = ("xt", "yt")
+    if src_crs is None:
+        src_crs = CRS.from_cf(ds.cf["grid_mapping"].attrs)
+    x, y = ds.cf["X"], ds.cf["Y"]
+    xt, yt = transform(x, y, src_crs, trg_crs)
+
+    return ds.assign_coords({trg_dims[0]: xt, trg_dims[1]: yt})
+
+
+def rotated_coord_transform(lon, lat, np_lon, np_lat, direction="rot2geo"):
+    """Transforms a coordinate into a rotated grid coordinate and vice versa.
+
+    The coordinates have to be given in degree and will be returned in degree.
+
+    Parameters
+    ----------
+    lon : float array like
+        Longitude coordinate.
+    lat : float array like
+        Latitude coordinate.
+    np_lon : float array like
+        Longitude coordinate of the rotated north pole.
+    np_lat : float array like
+        Latitude coordinate of the rotated north pole.
+    direction : str
+        Direction of the rotation.
+        Options are: 'rot2geo' (default) for a transformation to regular
+        coordinates from rotated. 'geo2rot' transforms regular coordinates
+        to rotated.
+
+    Returns
+    -------
+    lon_new : array like
+        New longitude coordinate.
+    lat_new : array like
+        New latitude coordinate.
+    """
+    warn(
+        "rotated_coord_transform is deprecated, please use transform_xy instead",
+        DeprecationWarning,
+        stacklevel=2,
+    )
+    # Convert degrees to radians
+    lon = np.deg2rad(lon)
+    lat = np.deg2rad(lat)
+
+    theta = 90.0 - np_lat  # Rotation around y-axis
+    phi = np_lon + 180.0  # Rotation around z-axis
+
+    # Convert degrees to radians
+    phi = np.deg2rad(phi)
+    theta = np.deg2rad(theta)
+
+    # Convert from spherical to cartesian coordinates
+    x = np.cos(lon) * np.cos(lat)
+    y = np.sin(lon) * np.cos(lat)
+    z = np.sin(lat)
+
+    # Regular -> Rotated
+    if direction == "geo2rot":
+        x_new = (
+            np.cos(theta) * np.cos(phi) * x
+            + np.cos(theta) * np.sin(phi) * y
+            + np.sin(theta) * z
+        )
+        y_new = -np.sin(phi) * x + np.cos(phi) * y
+        z_new = (
+            -np.sin(theta) * np.cos(phi) * x
+            - np.sin(theta) * np.sin(phi) * y
+            + np.cos(theta) * z
+        )
+
+    # Rotated -> Regular
+    elif direction == "rot2geo":
+        phi = -phi
+        theta = -theta
+
+        x_new = (
+            np.cos(theta) * np.cos(phi) * x
+            + np.sin(phi) * y
+            + np.sin(theta) * np.cos(phi) * z
+        )
+        y_new = (
+            -np.cos(theta) * np.sin(phi) * x
+            + np.cos(phi) * y
+            - np.sin(theta) * np.sin(phi) * z
+        )
+        z_new = -np.sin(theta) * x + np.cos(theta) * z
+
+    # Convert cartesian back to spherical coordinates
+    lon_new = np.arctan2(y_new, x_new)
+    lat_new = np.arcsin(z_new)
+
+    # Convert radians back to degrees
+    lon_new = np.rad2deg(lon_new)
+    lat_new = np.rad2deg(lat_new)
+
+    return lon_new, lat_new
+
+
+def grid_mapping(pollon, pollat, mapping_name=None):
+    """creates a grid mapping DataArray object"""
+    if mapping_name is None:
+        mapping_name = cf.DEFAULT_MAPPING_NCVAR
+    da = xr.DataArray(np.zeros((), dtype=np.int32))
+    attrs = cf.mapping.copy()
+    attrs["grid_north_pole_longitude"] = pollon
+    attrs["grid_north_pole_latitude"] = pollat
+    da.attrs = attrs
+    da.name = mapping_name
+    return da
diff --git a/cordex/core/utils.py b/cordex/core/utils.py
index 3e0e89a..ddf02bd 100644
--- a/cordex/core/utils.py
+++ b/cordex/core/utils.py
@@ -1,7 +1,4 @@
 import tempfile
-import warnings
-
-import xarray as xr
 
 from . import cf
 
@@ -28,106 +25,3 @@ def pole_crs(ds):
     from cartopy.crs import RotatedPole
 
     return RotatedPole(*pole(ds))
-
-
-def _map_crs(x_stack, y_stack, src_crs, trg_crs=None):
-    """coordinate transformation of longitude and latitude"""
-
-    from cartopy import crs as ccrs
-
-    if trg_crs is None:
-        trg_crs = ccrs.PlateCarree()
-    result = trg_crs.transform_points(src_crs, x_stack, y_stack)
-    return result[:, :, 0], result[:, :, 1]
-
-
-# wrapper function for xarray.apply_ufunc
-def map_crs(x, y, src_crs, trg_crs=None):
-    """coordinate transformation using cartopy
-
-    Transforms the coordinates x, y from the source crs
-    into the target crs using cartopy.
-
-    Parameters
-    ----------
-    x : float array like
-        x coordinate of source crs.
-    y : float array like
-        y coordinate of source crs.
-    src_crs : cartopy.crs
-        Source coordinate reference system in which x and y
-        are defined.
-    trg_crs : cartopy.crs
-        Target coordinate reference system into which x and y
-        should be transformed. If `None`, `PlateCarree` is used.
-
-    Returns
-    -------
-    x_map : xr.DataArray
-        Projected x coordinate.
-    y_map : xr.DataArray
-        Projected y coordinate.
-
-    """
-    warnings.warn("output shape has changed to apply to COARDS conventions")
-    y_stack, x_stack = xr.broadcast(y, x)
-    input_core_dims = 2 * [list(x_stack.dims)] + [[], []]
-    output_core_dims = 2 * [list(x_stack.dims)]
-    result = xr.apply_ufunc(
-        _map_crs,  # first the function
-        x_stack,  # now arguments in the order expected by 'interp1_np'
-        y_stack,
-        src_crs,
-        trg_crs,
-        input_core_dims=input_core_dims,  # list with one entry per arg
-        # [["rlat", "rlon"], ["rlat", "rlon"]],
-        output_core_dims=output_core_dims
-        # exclude_dims=set(("lat",)),  # dimensions allowed to change size. Must be set!
-    )
-    result[0].name = "x_map"
-    result[1].name = "y_map"
-    return result
-
-
-#
-#
-# def copy_dataset(src, varname=None, timestep=None, destination=None):
-#    """copy an existing NetCDF dataset on disk"""
-#    from netCDF4 import Dataset
-#
-#    if varname is None:
-#        variables = src.variables
-#    else:
-#        variables = {varname: src.variables[varname]}
-#    if destination is None:
-#        destination = "copy.nc"
-#    dst = Dataset(destination, "w")
-#    # copy attributes
-#    for name in src.ncattrs():
-#        dst.setncattr(name, getattr(src, name))
-#    # copy dimensions
-#    for name, dimension in src.dimensions.items():
-#        if timestep and name == "time":
-#            length = 1
-#        else:
-#            length = len(dimension) if not dimension.isunlimited() else None
-#        # dst.createDimension( name, (len(dimension) if not dimension.isunlimited() else None))
-#        dst.createDimension(name, length)
-#        # copy all file data except for the excluded
-#    for name, variable in variables.items():
-#        if hasattr(variable, "_FillValue"):
-#            fill_value = getattr(variable, "_FillValue")
-#        else:
-#            fill_value = None
-#        var = dst.createVariable(
-#            name, variable.dtype, variable.dimensions, fill_value=fill_value
-#        )
-#        for attr in variable.ncattrs():
-#            if attr != "_FillValue":
-#                dst.variables[name].setncattr(attr, getattr(variable, attr))
-#        if variable.shape:
-#            if timestep is None or "time" not in variable.dimensions:
-#                dst.variables[name][:] = src.variables[name][:]
-#            else:
-#                dst.variables[name][0] = src.variables[name][timestep]
-#    return dst
diff --git a/docs/api.rst b/docs/api.rst
index e9da5be..f5a3a5a 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -16,9 +16,9 @@ Top-level functions
    cordex_domain
    domain_info
    create_dataset
-   rotated_coord_transform
-   map_crs
    vertices
+   transform
+   transform_coords
 
 Tutorial
 ========
diff --git a/docs/whats_new.rst b/docs/whats_new.rst
index 7086cfd..050f2ed 100644
--- a/docs/whats_new.rst
+++ b/docs/whats_new.rst
@@ -16,10 +16,14 @@ New Features
 
 - Cmorization module from ``pyremo.cmor`` has been moved upstream into ``cordex.cmor``. A new cmorization funcion :py:meth:`cmor.cmorize_variable` now allows for easy cmorization using the archive specifications included in ``py-cordex`` (:pull:`68`, :pull:`74`, :pull:`76`, :pull:`77`).
 
+- Added :py:meth:`transform` in favour of deprecated :py:meth:`map_crs` and :py:meth:`rotated_coord_transform`. Dropped ``cartopy`` dependency in favour of ``pyproj`` (:pull:`71`).
+
 Internal Changes
 ~~~~~~~~~~~~~~~~
 
 - Updated CI pipeline (:pull:`83`).
+- Refactored transformation functions (:pull:`71`).
+- Pinned ``pyproj>=3.3.0`` and ``cf_xarray>=0.8.0`` (Dropped python 3.7 support).
 
 
 v0.4.1 (23 June 2022)


Test patch:
diff --git a/tests/test_domain.py b/tests/test_domain.py
index 3828a8c..b9d601c 100644
--- a/tests/test_domain.py
+++ b/tests/test_domain.py
@@ -1,13 +1,8 @@
-# -*- coding: utf-8 -*-
-# flake8: noqa
-
 import numpy as np
-import pytest
-import xarray as xr
 
 import cordex as cx
 
-from . import has_cartopy, requires_cartopy
+from . import requires_cartopy
 
 
 def test_domain_basic():
@@ -73,34 +68,7 @@ def test_domain_info():
     assert cx.domain_info("EUR-11", table) == info
 
 
-# @requires_cartopy
-def test_mapping():
-    import cartopy.crs as ccrs
-
-    eur11 = cx.cordex_domain("EUR-11")
-    pole = (
-        eur11.rotated_latitude_longitude.grid_north_pole_longitude,
-        eur11.rotated_latitude_longitude.grid_north_pole_latitude,
-    )
-    lon1, lat1 = cx.rotated_coord_transform(
-        eur11.rlon, eur11.rlat, *pole, direction="rot2geo"
-    )
-    transform = ccrs.RotatedPole(*pole)
-    lon2, lat2 = cx.map_crs(eur11.rlon, eur11.rlat, transform)
-
-    assert np.allclose(lon1.T, lon2)
-    assert np.allclose(lat1.T, lat2)
-
-    # test if retransforming of lon lat to rlon rlat gives correct results
-    rlon2, rlat2 = cx.map_crs(
-        eur11.lon, eur11.lat, src_crs=ccrs.PlateCarree(), trg_crs=transform
-    )
-    rlat1, rlon1 = xr.broadcast(eur11.rlat, eur11.rlon)
-
-    assert np.allclose(rlon1, rlon2)
-    assert np.allclose(rlat1, rlat2)
-
-
+@requires_cartopy
 def test_vertices():
     eur11 = cx.cordex_domain("EUR-11")
     import cartopy.crs as ccrs
@@ -110,3 +78,4 @@ def test_vertices():
         eur11.rotated_latitude_longitude.grid_north_pole_latitude,
     )
     vertices = cx.vertices(eur11.rlon, eur11.rlat, src_crs=ccrs.RotatedPole(*pole))
+    return vertices
diff --git a/tests/test_transform.py b/tests/test_transform.py
new file mode 100644
index 0000000..9b0a34f
--- /dev/null
+++ b/tests/test_transform.py
@@ -0,0 +1,57 @@
+import numpy as np
+import xarray as xr
+from cf_xarray.datasets import rotds
+from pyproj import CRS
+
+import cordex as cx
+
+from . import requires_cartopy
+
+
+@requires_cartopy
+def test_map_crs():
+    import cartopy.crs as ccrs
+
+    eur11 = cx.cordex_domain("EUR-11")
+    pole = (
+        eur11.rotated_latitude_longitude.grid_north_pole_longitude,
+        eur11.rotated_latitude_longitude.grid_north_pole_latitude,
+    )
+    lon1, lat1 = cx.rotated_coord_transform(
+        eur11.rlon, eur11.rlat, *pole, direction="rot2geo"
+    )
+    transform = ccrs.RotatedPole(*pole)
+    lon2, lat2 = cx.map_crs(eur11.rlon, eur11.rlat, transform)
+
+    assert np.allclose(lon1.T, lon2)
+    assert np.allclose(lat1.T, lat2)
+
+    # test if retransforming of lon lat to rlon rlat gives correct results
+    rlon2, rlat2 = cx.map_crs(
+        eur11.lon, eur11.lat, src_crs=ccrs.PlateCarree(), trg_crs=transform
+    )
+    rlat1, rlon1 = xr.broadcast(eur11.rlat, eur11.rlon)
+
+    assert np.allclose(rlon1, rlon2)
+    assert np.allclose(rlat1, rlat2)
+
+
+def test_transform():
+    ds = rotds.copy(deep=False)
+    rotated = CRS.from_cf(ds.cf["grid_mapping"].attrs)
+
+    xt, yt = cx.transform(ds.rlon, ds.rlat, rotated)
+    assert np.allclose(ds.lon, xt)
+    assert np.allclose(ds.lat, yt)
+
+    # test if retransforming of lon lat to rlon rlat gives correct results
+    rlon_t, rlat_t = cx.transform(xt, yt, CRS("EPSG:4326"), rotated)
+    rlat_b, rlon_b = xr.broadcast(ds.rlat, ds.rlon)
+
+    assert np.allclose(rlon_b, rlon_t)
+    assert np.allclose(rlat_b, rlat_t)
+
+    ds1 = cx.transform_coords(ds)
+
+    assert np.allclose(ds1.lon, ds1.xt)
+    assert np.allclose(ds1.lat, ds1.yt)