Add functionality for download of subregion (#70)

Co-authored-by: Aytaç PAÇAL <[email protected]>
Co-authored-by: Peter Kalverla <[email protected]>

README.rst

@@ -33,9 +33,10 @@ With era5cli you can: ​
 - select multiple variables for several months and years
 - split outputs by years, producing a separate file for every year instead of merging them in one file
 - download multiple files at once
+- extract data for a sub-region of the globe
 - download ERA5 back extension (preliminary version)
 
 .. inclusion-marker-end-do-not-remove
 
-| Free software: Apache Software License 2.0 
+| Free software: Apache Software License 2.0
 | Documentation: https://era5cli.readthedocs.io

era5cli/cli.py

@@ -159,6 +159,22 @@ def _build_parser():
                              ''')
     )
 
+    common.add_argument(
+        "--area", nargs=4, type=float,
+        required=False,
+        help=textwrap.dedent('''\
+                            Coordinates in case extraction of a subregion is
+                            requested. Specified as
+                            ymax xmin ymin xmax with x and y in the
+                            range -180, +180 and -90, +90, respectively
+                            e.g. --area 90 -180 -90 180.
+                            Requests are rounded down to two decimals.
+                            Without specification, the whole available area
+                            will be returned.
+
+                            ''')
+    )
+
     mnth = argparse.ArgumentParser(add_help=False)
 
     mnth.add_argument(
@@ -383,6 +399,7 @@ def _execute(args):
             days=days,
             hours=hours,
             variables=args.variables,
+            area=args.area,
             outputformat=args.format,
             outputprefix=args.outputprefix,
             period=args.command,

era5cli/fetch.py

@@ -26,6 +26,12 @@ class Fetch:
             hours to download data for.
         variables: list(str)
             List of variable names to download data for.
+        area: None, list(float)
+            Coordinates in case extraction of a subregion is requested.
+            Specified as [ymax, xmin, ymin, xmax], with x and y
+            in the range -180, +180 and -90, +90, respectively. Requests are
+            rounded down to two decimals. Without specification, the whole
+            available area will be returned.
         outputformat: str
             Type of file to download: 'netcdf' or 'grib'.
         outputprefix: str
@@ -65,7 +71,7 @@ class Fetch:
 
     def __init__(self, years: list, months: list, days: list,
                  hours: list, variables: list, outputformat: str,
-                 outputprefix: str, period: str, ensemble: bool,
+                 outputprefix: str, period: str, ensemble: bool, area=None,
                  statistics=None, synoptic=None, pressurelevels=None,
                  merge=False, threads=None, prelimbe=False, land=False):
         """Initialization of Fetch class."""
@@ -86,6 +92,9 @@ def __init__(self, years: list, months: list, days: list,
         """list(int): List of pressure levels."""
         self.variables = variables
         """list(str): List of variables."""
+        self.area = area
+        """list(float): Coordinates specifying the subregion that will be
+        extracted. Default is None for whole available area."""
         self.outputformat = outputformat
         """str: File format of output file."""
         self.years = years
@@ -149,6 +158,13 @@ def _extension(self):
             raise ValueError('Unknown outputformat: {}'.format(
                 self.outputformat))
 
+    def _process_areaname(self):
+        (ymax, xmin, ymin, xmax) = [round(c) for c in self.area]
+        def lon(x): return f"{x}E" if x >= 0 else f"{abs(x)}W"
+        def lat(y): return f"{y}N" if y >= 0 else f"{abs(y)}S"
+        name = f"_{lon(xmin)}-{lon(xmax)}_{lat(ymin)}-{lat(ymax)}"
+        return name
+
     def _define_outputfilename(self, var, years):
         """Define output filename."""
         start, end = years[0], years[-1]
@@ -157,6 +173,8 @@ def _define_outputfilename(self, var, years):
                   else self.outputprefix)
         yearblock = f"{start}-{end}" if not start == end else f"{start}"
         fname = f"{prefix}_{var}_{yearblock}_{self.period}"
+        if self.area:
+            fname += self._process_areaname()
         if self.ensemble:
             fname += "_ensemble"
         if self.statistics:
@@ -267,6 +285,30 @@ def _check_variable(self, variable):
                 "Invalid variable name: {}".format(variable)
             )
 
+    def _check_area(self):
+        """Confirm that area parameters are correct."""
+        (ymax, xmin, ymin, xmax) = self.area
+        if not (-90 <= ymax <= 90
+                and -90 <= ymin <= 90
+                and -180 <= xmin <= 180
+                and -180 <= xmax <= 180
+                and ymax > ymin
+                and xmax != xmin
+                ):
+            raise ValueError(
+                "Provide coordinates as ymax xmin ymin xmax. "
+                "x must be in range -180,+180 and y must be in range -90,+90."
+            )
+
+    def _parse_area(self):
+        """Parse area parameters to accepted coordinates."""
+        self._check_area()
+        area = [round(coord, ndigits=2) for coord in self.area]
+        if self.area != area:
+            print(
+                f"NB: coordinates {self.area} rounded down to two decimals.\n")
+        return area
+
     def _build_name(self, variable):
         """Build up name of dataset to use"""
 
@@ -310,6 +352,9 @@ def _build_request(self, variable, years):
             self._check_levels()
             request["pressure_level"] = self.pressure_levels
 
+        if self.area:
+            request["area"] = self._parse_area()
+
         product_type = self._product_type()
         if product_type is not None:
             request["product_type"] = product_type

tests/test_cli.py

@@ -28,6 +28,49 @@ def test_parse_args():
     assert not args.threads
     assert args.variables == ['total_precipitation']
     assert args.land
+    assert not args.area
+
+
+def test_area_argument():
+    """Test if area argument is parsed correctly."""
+    # Test if area arguments are parsed correctly
+    argv = ['hourly', '--startyear', '2008',
+            '--variables', 'total_precipitation', '--statistics',
+            '--endyear', '2008', '--ensemble',
+            '--area', '90', '-180', '-90', '180']
+    args = cli._parse_args(argv)
+    assert args.area == [90, -180, -90, 180]
+
+    # Check that area defaults to None
+    argv = ['hourly', '--startyear', '2008',
+            '--variables', 'total_precipitation', '--statistics',
+            '--endyear', '2008', '--ensemble']
+    args = cli._parse_args(argv)
+    assert not args.area
+
+    # Requires four values
+    with pytest.raises(SystemExit):
+        argv = ['hourly', '--startyear', '2008',
+                '--variables', 'total_precipitation', '--statistics',
+                '--endyear', '2008', '--ensemble',
+                '--area', '90', '-180', '-90']
+        cli._parse_args(argv)
+
+    # A value cannot be missing
+    with pytest.raises(SystemExit):
+        argv = ['hourly', '--startyear', '2008',
+                '--variables', 'total_precipitation', '--statistics',
+                '--endyear', '2008', '--ensemble',
+                '--area', '90', '-180', '-90', '']
+        cli._parse_args(argv)
+
+    # Values must be numeric
+    with pytest.raises(SystemExit):
+        argv = ['hourly', '--startyear', '2008',
+                '--variables', 'total_precipitation', '--statistics',
+                '--endyear', '2008', '--ensemble',
+                '--area', '90', '-180', '-90', 'E']
+        cli._parse_args(argv)
 
 
 def test_period_args():

tests/test_fetch.py

@@ -7,7 +7,8 @@
 
 def initialize(outputformat='netcdf', merge=False, statistics=None,
                synoptic=None, ensemble=True, pressurelevels=None,
-               threads=2, period='hourly', variables=['total_precipitation'],
+               threads=2, period='hourly', area=None,
+               variables=['total_precipitation'],
                years=[2008, 2009], months=list(range(1, 13)),
                days=list(range(1, 32)), hours=list(range(0, 24)),
                prelimbe=False, land=False):
@@ -16,6 +17,7 @@ def initialize(outputformat='netcdf', merge=False, statistics=None,
                        months=months,
                        days=days,
                        hours=hours,
+                       area=area,
                        variables=variables,
                        outputformat=outputformat,
                        outputprefix='era5',
@@ -201,6 +203,27 @@ def test_define_outputfilename():
     fn = 'era5-land_total_precipitation_2008_hourly.nc'
     assert fname == fn
 
+    era5.area = [90.0, -180.0, -90.0, 180.0]
+    fname = era5._define_outputfilename('total_precipitation', [2008])
+    fn = (
+        'era5-land_total_precipitation_2008_hourly_180W-180E_90S-90N.nc'
+    )
+    assert fname == fn
+
+    era5.area = [90.0, -180.0, -80.999, 170.001]
+    fname = era5._define_outputfilename('total_precipitation', [2008])
+    fn = (
+        'era5-land_total_precipitation_2008_hourly_180W-170E_81S-90N.nc'
+    )
+    assert fname == fn
+
+    era5.area = [0, 120, -90, 180]
+    fname = era5._define_outputfilename('total_precipitation', [2008])
+    fn = (
+        'era5-land_total_precipitation_2008_hourly_120E-180E_90S-0N.nc'
+    )
+    assert fname == fn
+
 
 def test_number_outputfiles(capsys):
     """Test function for the number of outputs."""
@@ -520,3 +543,58 @@ def test_more_incompatible_options():
     era5 = initialize(statistics=True, ensemble=False)
     with pytest.raises(ValueError):
         era5._build_request('total_precipitation', [2008])
+
+
+def test_area():
+    """Test that area is parsed properly."""
+    era5 = initialize()
+    assert era5.area is None
+
+    era5 = initialize(area=[90, -180, -90, 180])
+    (name, request) = era5._build_request('total_precipitation', [2008])
+    assert era5.area == [90, -180, -90, 180]
+    assert request["area"] == [90, -180, -90, 180]
+
+    # Decimals are rounded down
+    era5 = initialize(area=[89.9999, -179.90, -90.0000, 179.012])
+    (name, request) = era5._build_request('total_precipitation', [2008])
+    assert request["area"] == [90.0, -179.90, -90.0, 179.01]
+
+    # ymax may not be lower than ymin
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[-10, -180, 10, 180])
+        era5._build_request('total_precipitation', [2008])
+
+    # xmin higher than xmax should be ok
+    era5 = initialize(area=[90, 120, -90, -120])
+    (name, request) = era5._build_request('total_precipitation', [2008])
+    assert request["area"] == [90.0, 120.0, -90.0, -120.0]
+
+    # ymax may not equal ymin
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[0, -180, 0, 180])
+        era5._build_request('total_precipitation', [2008])
+
+    # xmin may not equal xmax
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[90, 0, -90, 0])
+        era5._build_request('total_precipitation', [2008])
+
+    # ymax, xmin, ymin, xmax may not be out of bounds
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[1000, -180, -90, 180])
+        era5._build_request('total_precipitation', [2008])
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[90, 1000, -90, 180])
+        era5._build_request('total_precipitation', [2008])
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[90, -180, 1000, 180])
+        era5._build_request('total_precipitation', [2008])
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[90, -180, -90, 1000])
+        era5._build_request('total_precipitation', [2008])
+
+    # Coordinate missing
+    with pytest.raises(ValueError):
+        era5 = initialize(area=[-180, 180, -90])
+        era5._build_request('total_precipitation', [2008])