added function to add 2D latlon to dataset

.travis.yml

@@ -31,7 +31,7 @@ before_install:
 - conda config --add channels conda-forge
 - conda config --set channel_priority strict
 # Create conda environment
-- conda create -n test python=$PYTHON_VERSION rasterio=1.0.* scipy xarray netcdf4 dask
+- conda create -n test python=$PYTHON_VERSION rasterio=1.0.* scipy xarray pyproj netcdf4 dask
 - source activate test
 
 install:

appveyor.yml

@@ -16,7 +16,7 @@ install:
 #-----------------------------------------------------------------------------
 # Create conda environment
 #-----------------------------------------------------------------------------
-  - conda create -n test python=%PYTHON_VERSION% rasterio=1.0.* scipy xarray netcdf4 dask
+  - conda create -n test python=%PYTHON_VERSION% rasterio=1.0.* scipy xarray pyproj netcdf4 dask
   - activate test
 
 #-------------------------------------------------------------------------------

rioxarray/rioxarray.py

@@ -13,8 +13,11 @@
 import copy
 from abc import abstractmethod
 from datetime import datetime
+from distutils.version import LooseVersion
+from functools import partial
 
 import numpy as np
+import pyproj
 import rasterio.warp
 import xarray
 from affine import Affine
@@ -36,6 +39,7 @@
 FILL_VALUE_NAMES = ("_FillValue", "missing_value", "fill_value", "nodata")
 UNWANTED_RIO_ATTRS = ("nodatavals", "crs", "is_tiled", "res")
 DEFAULT_GRID_MAP = "spatial_ref"
+_PYPROJ22 = LooseVersion(pyproj.__version__) >= LooseVersion("2.2.0")
 
 
 def affine_to_coords(affine, width, height, x_dim="x", y_dim="y"):
@@ -200,6 +204,32 @@ def _make_dst_affine(src_data_array, src_crs, dst_crs, dst_resolution=None):
     return dst_affine, dst_width, dst_height
 
 
+def _get_proj_transformer(src_crs, dst_crs):
+    """
+    Get pyproj Transformer based on version of pyproj
+
+    Parameters
+    ----------
+    src_crs: input to create source pyproj.CRS/pyproj.Proj
+    src_crs: input to create destination pyproj.CRS/pyproj.Proj
+
+    Returns
+    -------
+    pyproj Transformer
+    """
+    if _PYPROJ22:
+        return pyproj.Transformer.from_crs(src_crs, dst_crs, always_xy=True).transform
+
+    def ensure_proj(in_crs):
+        if not isinstance(in_crs, pyproj.Proj):
+            if hasattr(in_crs, "to_dict"):
+                in_crs = in_crs.to_dict()
+            in_crs = pyproj.Proj(in_crs, preserve_units=True)
+        return in_crs
+
+    return partial(pyproj.transform, ensure_proj(src_crs), ensure_proj(dst_crs))
+
+
 class XRasterBase(object):
     """This is the base class for the GIS extensions for xarray"""
 
@@ -390,6 +420,54 @@ def shape(self):
         """tuple: Returns the shape (width, height)"""
         return (self.width, self.height)
 
+    def add_latlon(self, inplace=False):
+        """Adds 2D 'latitude' and 'longitude' coordinates to
+        an `xarray.Dataset` or `xarray.DataArray` based on the
+        projection and the 'x' and 'y' coordinates. If the 'latitude'
+        and 'longitude' coordinates already exist, it will do nothing.
+
+        Parameters
+        ----------
+        inplace: bool
+            If True, it will perform the operation inplace. Default is False.
+
+        Returns
+        -------
+        :obj:`xarray.DataArray` or :obj:`xarray.Dataset`
+            The data_array with the dimensions added.
+        """
+        # don't do anything if the coordinates exist already
+        if all(coord in self._obj.coords for coord in ("latitude", "longitude")):
+            return self._obj
+
+        # Compute the lon/lat coordinates with pyproj.transform
+        transformer = _get_proj_transformer(self.crs, CRS.from_epsg(4326))
+        x_coords, y_coords = np.meshgrid(
+            self._obj.coords[self.x_dim], self._obj.coords[self.y_dim]
+        )
+        lon, lat = transformer(x_coords.flatten(), y_coords.flatten())
+        lon = np.asarray(lon).reshape((self.height, self.width))
+        lat = np.asarray(lat).reshape((self.height, self.width))
+
+        if inplace:
+            xds = self._obj
+        else:
+            xds = self._obj.copy()
+        # add coordinates to dataset
+        xds.coords["longitude"] = ((self.y_dim, self.x_dim), lon)
+        xds.coords["longitude"].attrs = {
+            "long_name": "longitude",
+            "standard_name": "longitude",
+            "units": "degrees_east",
+        }
+        xds.coords["latitude"] = ((self.y_dim, self.x_dim), lat)
+        xds.coords["latitude"].attrs = {
+            "long_name": "latitude",
+            "standard_name": "latitude",
+            "units": "degrees_north",
+        }
+        return xds
+
 
 @xarray.register_dataarray_accessor("rio")
 class RasterArray(XRasterBase):

setup.py

@@ -24,7 +24,7 @@ def get_version():
 with open("README.rst") as readme_file:
     readme = readme_file.read()
 
-requirements = ["rasterio", "scipy", "xarray"]
+requirements = ["rasterio", "scipy", "xarray", "pyproj"]
 
 test_requirements = ["pytest>=3.6", "pytest-cov", "mock"]
 

sphinx/history.rst

@@ -4,6 +4,7 @@ History
 0.0.10
 -----
 - Add support for opening netcdf/hdf files with `rioxarray.open_rasterio` (issue #32)
+- Added `add_latlon()` method for 2D latlon coordinates (pull #34)
 
 0.0.9
 -----

test/integration/test_integration_rioxarray.py

@@ -975,3 +975,70 @@ def test_crs_writer__missing():
         test_da.rio.write_crs()
     with pytest.raises(MissingCRS):
         test_da.to_dataset(name="test").rio.write_crs()
+
+
+@pytest.mark.parametrize(
+    "input_ds,compare_ds,xarray_open",
+    [
+        (
+            os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"),
+            os.path.join(TEST_COMPARE_DATA_DIR, "MODIS_ARRAY_LATLON.nc"),
+            xarray.open_dataarray,
+        ),
+        (
+            os.path.join(TEST_INPUT_DATA_DIR, "MODIS_MYD09GQ_15TWF79_20171019.nc"),
+            os.path.join(
+                TEST_COMPARE_DATA_DIR, "MODIS_MYD09GQ_15TWF79_20171019_LATLON.nc"
+            ),
+            xarray.open_dataset,
+        ),
+    ],
+)
+def test_latlon(input_ds, compare_ds, xarray_open):
+    with xarray_open(input_ds, autoclose=True) as xds, xarray_open(
+        compare_ds, autoclose=True
+    ) as xdc:
+        xds_output = xds.rio.add_latlon()
+        assert "latitude" not in xds.coords
+        assert "longitude" not in xds.coords
+        xarray.testing.assert_allclose(xds_output, xdc)
+
+
+@pytest.mark.parametrize(
+    "input_ds,compare_ds,xarray_open",
+    [
+        (
+            os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"),
+            os.path.join(TEST_COMPARE_DATA_DIR, "MODIS_ARRAY_LATLON.nc"),
+            xarray.open_dataarray,
+        ),
+        (
+            os.path.join(TEST_COMPARE_DATA_DIR, "MODIS_ARRAY_LATLON.nc"),
+            os.path.join(TEST_COMPARE_DATA_DIR, "MODIS_ARRAY_LATLON.nc"),
+            xarray.open_dataarray,
+        ),
+        (
+            os.path.join(TEST_INPUT_DATA_DIR, "MODIS_MYD09GQ_15TWF79_20171019.nc"),
+            os.path.join(
+                TEST_COMPARE_DATA_DIR, "MODIS_MYD09GQ_15TWF79_20171019_LATLON.nc"
+            ),
+            xarray.open_dataset,
+        ),
+        (
+            os.path.join(
+                TEST_COMPARE_DATA_DIR, "MODIS_MYD09GQ_15TWF79_20171019_LATLON.nc"
+            ),
+            os.path.join(
+                TEST_COMPARE_DATA_DIR, "MODIS_MYD09GQ_15TWF79_20171019_LATLON.nc"
+            ),
+            xarray.open_dataset,
+        ),
+    ],
+)
+def test_latlon__inplace(input_ds, compare_ds, xarray_open):
+    with xarray_open(input_ds, autoclose=True) as xds, xarray_open(
+        compare_ds, autoclose=True
+    ) as xdc:
+        out_xds = xds.rio.add_latlon(inplace=True)
+        xarray.testing.assert_allclose(out_xds, xdc)
+        xarray.testing.assert_allclose(xds, xdc)