Implement smecv-grid (#6)

AUTHORS.rst

@@ -2,5 +2,5 @@
 Developers
 ==========
 
-* tscanlon <[email protected]>
-* wpreimes <[email protected]>
+* Tracy Scanlon <[email protected]>
+* Wolfgang Preimesberger <[email protected]>

CHANGES.rst

@@ -4,7 +4,12 @@ Changelog
 
 Version 0.X (unreleased)
 ========================
+-
 
+Version 0.1
+========================
+
+- Add smecv grid support
 - Add reader from images and image stacks
 - Add reshuffle module
 - Add first version of metadata module

c3s_sm/grid.py

@@ -1,71 +1,11 @@
 # -*- coding: utf-8 -*-
 
-from pygeogrids.grids import BasicGrid, CellGrid
-import numpy as np
-try:
-    from smecv_grid.grid import SMECV_Grid_v042
-except:
-    import warnings
-    warnings.warn('SMECV grid is not installed')
+from smecv_grid.grid import SMECV_Grid_v042
 
 def C3SCellGrid():
-    '''
-    Returns
-    -------
-    grid : CellGrid
-        The global QDEG grid
-    '''
-    resolution = 0.25
-    lon, lat = np.meshgrid(
-        np.arange(-180 + resolution / 2, 180 + resolution / 2, resolution),
-        np.flipud(np.arange(-90 + resolution / 2, 90 + resolution / 2, resolution)))
-
-    return BasicGrid(lon.flatten(), lat.flatten()).to_cell_grid(cellsize=5.)
-
-def C3SLandPoints(grid):
-    '''
-    Create a subset of land points.
-    Returns
-    -------
-    points : np.array
-        Points in the passed grid that are over land
-    dist : np.array
-        Distance between the reference land points and the next point in the
-        passed grid.
-    '''
-    lg = SMECV_Grid_v042('land')
-    lat = lg.get_grid_points()[2]
-    lon = lg.get_grid_points()[1]
-
-    points, dist =  grid.find_nearest_gpi(lon, lat)
-
-    return points, dist
-
+    return SMECV_Grid_v042(None)
 
 def C3SLandGrid():
-    '''
-    0.25deg cell grid of land points from c3s land mask.
-
-    Returns
-    -------
-    landgrid : CellGrid
-        The reduced QDEG grid
-    '''
-    grid = C3SCellGrid()
-    land_gpis, dist = C3SLandPoints(grid)
-    if any(dist) > 0:
-        raise Exception('C3S grid does not conform with QDEG grid')
-    return grid.subgrid_from_gpis(land_gpis)
-
-
-if __name__ == '__main__':
-    grid = C3SCellGrid()
-    shape = grid.get_grid_points()[0].size
-    assert shape == 1036800
-
-    landpoints, dist = C3SLandPoints(grid)
+    return SMECV_Grid_v042('land')
 
-    land_grid = C3SLandGrid()
-    shape = land_grid.get_grid_points()[0].size
-    assert shape == 244243
 

c3s_sm/interface.py

@@ -188,9 +188,9 @@ def read(self, timestamp=None):
             for attr in variable.ncattrs():
                 param_metadata.update({str(attr): getattr(variable, attr)})
 
-            param_data = variable[0][:].filled().flatten()
+            param_data = np.flipud(variable[0][:].filled()).flatten()
 
-            param_img[str(param)] = param_data[np.sort(self.grid.activegpis)]
+            param_img[str(param)] = param_data[self.grid.activegpis]
             img_meta[param] = param_metadata
 
         # add global attributes
@@ -203,11 +203,12 @@ def read(self, timestamp=None):
             return Image(self.grid.activearrlon, self.grid.activearrlat,
                          param_img, img_meta, timestamp)
         else:
+            yres, xres = self.grid.shape
             for key in param_img:
-                param_img[key] = param_img[key].reshape(720, 1440)
+                param_img[key] = param_img[key].reshape(xres, yres)
 
-            return Image(self.grid.activearrlon.reshape(720, 1440),
-                         self.grid.activearrlat.reshape(720, 1440),
+            return Image(self.grid.activearrlon.reshape(xres, yres),
+                         self.grid.activearrlat.reshape(xres, yres),
                          param_img,
                          img_meta,
                          timestamp)

docs/img2ts.rst

@@ -13,7 +13,7 @@ methods. We have chosen to do it in the following way:
   `Orthogonal multidimensional array representation
   <http://cfconventions.org/cf-conventions/v1.6.0/cf-conventions.html#_orthogonal_multidimensional_array_representation>`_
 - Store the time series in 5x5 degree cells. This means there will be 2566 cell
-  files (without reduction to land points) and a file called ``grid.nc``
+  files (1001 when reduced to land points) and a file called ``grid.nc``
   which contains the information about which grid point is stored in which file.
   This allows us to read a whole 5x5 degree area into memory and iterate over the time series quickly.
 
@@ -25,10 +25,11 @@ program. An example would be:
 
 .. code-block:: shell
 
-   c3s_repurpose /c3s_images /timeseries/data 2000-01-01 2001-01-01 sm sm_uncertainty
+   c3s_repurpose /c3s_images /timeseries/data 2000-01-01 2001-01-01 sm sm_uncertainty --land_points True
 
 Which would take C3S SM data stored in ``/c3s_images`` from January 1st
-2000 to January 1st 2001 and store the parameters for soil moisture and its uncertainty as time
+2000 to January 1st 2001 and store the parameters for soil moisture and its uncertainty
+of points marked as 'land' in the smecv-grid as time
 series in the folder ``/timeseries/data``.
 
 Conversion to time series is performed by the `repurpose package

docs/index.rst

@@ -5,6 +5,8 @@
 
 .. include:: img2ts.rst
 
+.. include:: varnames.rst
+
 Contents
 ========
 

docs/varnames.rst

@@ -18,7 +18,7 @@ C3S SM variables as in the netcdf image files (and time series from netcdf image
 +-------------+---------------------------------+-------------------------------------------+
 | sensor      | Sensor Flag                     |                                           |
 +-------------+---------------------------------+-------------------------------------------+
-| sm          | Volumetric Soil Moisture        | [m3 m-3"]                                 |
+| sm          | Volumetric Soil Moisture        | [m3 m-3]                                 |
 +-------------+---------------------------------+-------------------------------------------+
-| time        | Time                            | ["days since 1970-01-01 00:00:00 UTC"]    |
+| time        | Time                            | [days since 1970-01-01 00:00:00 UTC]    |
 +-------------+---------------------------------+-------------------------------------------+

environment.yml

@@ -5,14 +5,13 @@ dependencies:
   - netcdf4=1.2.2
   - numpy
   - pandas
-  - scipy
   - pytest=3.1.2
   - pip:
     - pynetcf
     - repurpose
-    - pytesmo
     - pyresample
     - parse
     - pytest-cov==2.2.1
-    # - smecv_grid #### when public
+    - smecv_grid
+
 

requirements.txt

@@ -2,9 +2,12 @@
 # numpy
 # scipy>=0.9
 numpy
+pandas
 pygeobase
 netcdf4==1.2.2
 pyresample
 repurpose
 pynetcf
+repurpose
 parse
+smecv_grid

setup.cfg

@@ -1,10 +1,10 @@
 [metadata]
 name = c3s_sm
-summary = Readers and coverters for the C3S Soil Moisture data set
-author = tscanlon
-author-email = tracy.scanlon@geo.tuwien.ac.at
-license = none
-home-page = http://...
+summary = Readers and time series coverters for the C3S Soil Moisture data set
+author = Wolfgang Preimesberger
+author-email = wolfgang.preimesberger@geo.tuwien.ac.at
+license = MIT
+home-page = https://climate.copernicus.eu/
 description-file = README.rst
 # Add here all kinds of additional classifiers as defined under
 # https://pypi.python.org/pypi?%3Aaction=list_classifiers
@@ -15,7 +15,7 @@ classifier =
 [entry_points]
 # Add here console scripts like:
 console_scripts =
-     c3s_reshuffle = c3s_sm.reshuffle:run
+     c3s_repurpose = c3s_sm.reshuffle:run
 # For example:
 # console_scripts =
 #     fibonacci = esa_cci_sm.skeleton:run

tests/conftest.py

@@ -1,12 +1,11 @@
 #!/usr/bin/env python
 # -*- coding: utf-8 -*-
 """
-    Dummy conftest.py for c3s_sm_reader.
+    Dummy conftest.py for c3s_sm.
 
     If you don't know what this is for, just leave it empty.
     Read more about conftest.py under:
     https://pytest.org/latest/plugins.html
 """
 from __future__ import print_function, absolute_import, division
-
 import pytest

tests/test_grid.py

@@ -1,10 +1,39 @@
 # -*- coding: utf-8 -*-
 
-def test_cell_grid():
-    pass
+from c3s_sm.grid import C3SCellGrid
+from c3s_sm.grid import C3SLandGrid
+import numpy as np
 
 
-def test_landgrid():
-    pass
+def test_C3SCellGrid():
+    grid = C3SCellGrid()
+    gp, dist = grid.find_nearest_gpi(75.625, 14.625)
+    assert gp == 602942
+    lon, lat = grid.gpi2lonlat(602942)
+    assert lon == 75.625
+    assert lat == 14.625
+    assert np.where(grid.get_grid_points()[0] == 602942)[0][0] == 434462 # index
+    assert grid.get_grid_points()[1][434462] == lon
+    assert grid.get_grid_points()[2][434462] == lat
+    assert grid.gpi2cell(602942) == 1856
+    assert grid.gpis.size == 1036800
+    assert grid.gpis[0] == 1035360
+    assert np.unique(grid.get_grid_points()[3]).size == 2592
 
 
+def test_landgrid():
+    grid = C3SLandGrid()
+    gp, dist = grid.find_nearest_gpi(75.625, 14.625)
+    assert gp == 602942
+    lon, lat = grid.gpi2lonlat(602942)
+    assert lon == 75.625
+    assert lat == 14.625
+    assert np.where(grid.get_grid_points()[0] == 602942)[0][0] == 177048  # index
+    assert grid.get_grid_points()[1][177048] == lon
+    assert grid.get_grid_points()[2][177048] == lat
+    assert grid.gpi2cell(602942) == 1856
+    assert grid.gpis.size == 1036800
+    assert grid.activegpis.size == 244243
+    assert grid.gpis[0] == 1035360
+    assert grid.activegpis[0] == 999942
+    assert np.unique(grid.get_grid_points()[3]).size == 1001

tests/test_interface_img.py

@@ -3,6 +3,7 @@
 from c3s_sm.interface import C3SImg
 import os
 import numpy.testing as nptest
+from c3s_sm.grid import C3SLandGrid, C3SCellGrid
 
 # lat=48.125, lon=16.375
 def test_C33Ts_tcdr_combined_daily():
@@ -96,7 +97,51 @@ def test_C33Ts_icdr_passive_decadal():
     #assert(image.metadata['creator_name'] == 'Earth Observation Data Center (EODC)')
 
 
+def test_1Dreading():
+    ''' Test 1D reading with and without land grid, and if the results are the same'''
+    file = os.path.join(os.path.join(os.path.dirname(__file__),
+                        'c3s_sm-test-data', 'img', 'ICDR', '060_dailyImages', 'combined', '2017',
+                        'C3S-SOILMOISTURE-L3S-SSMV-COMBINED-DAILY-20170701000000-ICDR-v201706.0.0.nc'))
+
+    ds = C3SImg(file, mode='r', parameters='sm', array_1D=True)
+    image = ds.read()
+
+    assert ds.grid.find_nearest_gpi(75.625, 14.625) == (602942, 0)
+    ref_sm = image.data['sm'][434462]
+    ref_lat = image.lat[434462]
+    ref_lon = image.lon[434462]
+
+    assert ref_lat == 14.625
+    assert ref_lon == 75.625
+    nptest.assert_almost_equal(ref_sm, 0.360762, 5)
+    assert(image.metadata['sm']['long_name'] == 'Volumetric Soil Moisture')
+
+    land_grid = C3SLandGrid()
+
+    ds = C3SImg(file, mode='r', parameters='sm', array_1D=True, subgrid=land_grid)
+    image = ds.read()
+
+    assert ds.grid.find_nearest_gpi(75.625, 14.625) == (602942, 0)
+
+    sm = image.data['sm'][177048]
+    lat = image.lat[177048]
+    lon = image.lon[177048]
+
+    assert ref_lat == lat
+    assert ref_lon == lon
+    nptest.assert_almost_equal(ref_sm, sm, 5)
+
+    assert(image.metadata['sm']['long_name'] == 'Volumetric Soil Moisture')
+
+
+
+
+
+
+
+
 if __name__ == '__main__':
+    test_1Dreading()
     test_C33Ts_tcdr_combined_daily()
     test_C33Ts_tcdr_active_monthly()
     test_C33Ts_tcdr_passive_decadal()

tests/test_reshuffle.py

@@ -5,7 +5,6 @@
 import tempfile
 import numpy as np
 import numpy.testing as nptest
-import shutil
 
 from c3s_sm.reshuffle import main, parse_filename
 from c3s_sm.interface import C3STs
@@ -37,14 +36,14 @@ def test_reshuffle_TCDR_daily():
     enddate = "1991-08-08"
     #parameters = ["sm", "sm_uncertainty", "dnflag", "flag", "freqbandID", "mode", "sensor", "t0"]
     parameters = ['--parameters', 'sm', 'sm_uncertainty']
-    land_points = 'False' # TODO: set this to True when the smecv_grid is public
+    land_points = 'True'
 
     ts_path = tempfile.mkdtemp()
     args = [inpath, ts_path, startdate, enddate]  + \
            parameters + ['--land_points', land_points]
     main(args)
 
-    assert len(glob.glob(os.path.join(ts_path, "*.nc"))) == 2593 # todo: set this to 1002 when using land points only
+    assert len(glob.glob(os.path.join(ts_path, "*.nc"))) == 1002
 
     ds = C3STs(ts_path, remove_nans=True)
     ts = ds.read(75.625, 14.625)
@@ -55,6 +54,8 @@ def test_reshuffle_TCDR_daily():
                                        dtype=np.float32)
     nptest.assert_allclose(ts['sm_uncertainty'].values, ts_uncert_values_should,rtol=1e-5)
 
+    nptest.assert_almost_equal(ts['sm'].values, ds.read(602942)['sm'].values)
+
     #nptest.assert_allclose(ts['dnflag'].values, [1,0,2,1],rtol=1e-5)
     #nptest.assert_allclose(ts['flag'].values, [0,127,0,0],rtol=1e-5)
     #nptest.assert_allclose(ts['freqbandID'].values, [2,0,2,2],rtol=1e-5)