diff --git a/AUTHORS.rst b/AUTHORS.rst index 15eb050..cdc92e5 100644 --- a/AUTHORS.rst +++ b/AUTHORS.rst @@ -2,5 +2,5 @@ Developers ========== -* tscanlon -* wpreimes +* Tracy Scanlon +* Wolfgang Preimesberger diff --git a/CHANGES.rst b/CHANGES.rst index b1f8c48..4e28e73 100644 --- a/CHANGES.rst +++ b/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 diff --git a/c3s_sm/grid.py b/c3s_sm/grid.py index 73015eb..5463c97 100644 --- a/c3s_sm/grid.py +++ b/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 diff --git a/c3s_sm/interface.py b/c3s_sm/interface.py index 57c10d6..ffaf7f7 100644 --- a/c3s_sm/interface.py +++ b/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) diff --git a/docs/img2ts.rst b/docs/img2ts.rst index 316965d..5a9821a 100644 --- a/docs/img2ts.rst +++ b/docs/img2ts.rst @@ -13,7 +13,7 @@ methods. We have chosen to do it in the following way: `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 diff --git a/docs/index.rst b/docs/index.rst index c910f2b..3c130dc 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -5,6 +5,8 @@ .. include:: img2ts.rst +.. include:: varnames.rst + Contents ======== diff --git a/docs/varnames.rst b/docs/varnames.rst index 5100262..87a3a72 100644 --- a/docs/varnames.rst +++ b/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] | +-------------+---------------------------------+-------------------------------------------+ \ No newline at end of file diff --git a/environment.yml b/environment.yml index 39bf5ef..0acc6f5 100644 --- a/environment.yml +++ b/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 + diff --git a/requirements.txt b/requirements.txt index c9a33c9..6d592c8 100644 --- a/requirements.txt +++ b/requirements.txt @@ -2,9 +2,12 @@ # numpy # scipy>=0.9 numpy +pandas pygeobase netcdf4==1.2.2 pyresample repurpose pynetcf +repurpose parse +smecv_grid \ No newline at end of file diff --git a/setup.cfg b/setup.cfg index 1eb399b..7d9b544 100644 --- a/setup.cfg +++ b/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 diff --git a/tests/conftest.py b/tests/conftest.py index 5726edf..a0d834b 100644 --- a/tests/conftest.py +++ b/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 diff --git a/tests/test_grid.py b/tests/test_grid.py index 02788cd..111fe9c 100644 --- a/tests/test_grid.py +++ b/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 diff --git a/tests/test_interface_img.py b/tests/test_interface_img.py index ff64027..762f6c7 100644 --- a/tests/test_interface_img.py +++ b/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() diff --git a/tests/test_reshuffle.py b/tests/test_reshuffle.py index 7fc983b..ef7bdbc 100644 --- a/tests/test_reshuffle.py +++ b/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)