Adding the command line interface for shape subsetting and documentin…

…g the option

cis/parse.py

@@ -522,7 +522,10 @@ def get_subset_limits(subsetlimits, parser):
         # or
         # <dim_name>=[<start_value>]
         match = re.match(r'(?P<dim>[^=]+)=\[(?P<start>[^],]+)(?:,(?P<end>[^],]+))?\]$', seg)
-        if match is None or match.group('dim') is None or match.group('start') is None:
+        if seg.startswith('shape'):
+            # Don't use the regexp for this as it gets confused with the commas
+            limit_dict['shape'] = seg.split('=')[1]
+        elif match is None or match.group('dim') is None or match.group('start') is None:
             parser.error(
                 "A dimension for subsetting does not have dimension name, start value and/or end value specified")
         else:
@@ -629,7 +632,7 @@ def parse_colon_and_comma_separated_arguments(inputs, parser, options, compulsor
     return input_dicts
 
 
-def split_outside_brackets(input, seps=[','], brackets={'[': ']'}):
+def split_outside_brackets(input, seps=[','], brackets={'[': ']', '(': ')'}):
     """Splits an input string at separators that are not within brackets.
     :param input: input string to parse
     :param seps: list of separator characters - default: comma

cis/subsetting/subset.py

@@ -288,7 +288,8 @@ def _get_indices_for_lat_lon_points(lats, lons, region):
 
 
 def _get_ungridded_subset_region_indices(ungridded_data, region):
-    return _get_indices_for_lat_lon_points(ungridded_data.lon.data.flat, ungridded_data.lat.data.flat, region)
+    # We have to use flatten rather than flat, GEOS creates a copy of the data if it's a view anyway.
+    return _get_indices_for_lat_lon_points(ungridded_data.lon.data.flatten(), ungridded_data.lat.data.flatten(), region)
 
 
 def _get_gridded_subset_region_indices(gridded_data, region):

cis/test/integration/test_subset.py

@@ -25,6 +25,19 @@ def test_GIVEN_single_variable_in_ungridded_file_WHEN_subset_THEN_subsetted_corr
         self.check_latlon_subsetting(lat_max, lat_min, lon_max, lon_min)
         self.check_output_contains_variables(self.OUTPUT_FILENAME, [variable])
 
+    def test_GIVEN_shape_WHEN_subset_ungridded_data_THEN_subsetted_correctly(self):
+        variable = valid_aerosol_cci_variable
+        filename = valid_aerosol_cci_filename
+        lon_min, lon_max = -10, 10
+        lat_min, lat_max = 40, 60
+        shape_wkt = "POLYGON((-10 50, 0 60, 10 50, 0 40, -10 50))"
+        arguments = ['subset', variable + ':' + escape_colons(filename),
+                     'shape=%s' % shape_wkt, '-o', self.OUTPUT_FILENAME]
+        main_arguments = parse_args(arguments)
+        subset_cmd(main_arguments)
+        self.check_latlon_subsetting(lat_max, lat_min, lon_max, lon_min)
+        self.check_output_contains_variables(self.OUTPUT_FILENAME, [variable])
+
     def test_GIVEN_single_variable_as_var_name_in_ungridded_file_WHEN_subset_THEN_subsetted_correctly(self):
         variable = valid_aerosol_cci_variable
         filename = valid_aerosol_cci_filename

doc/subsetting.rst

@@ -41,6 +41,11 @@ where:
       using ``x=[90,270]``. The range between the start and end must not be greater than 360 degrees. The output
       coordinates will be on the requested grid, not the grid of the source data.
 
+    .. note::
+      An arbitrary lat/lon shape can also be provided using the ``shape`` limit and passing a valid WKT string as the
+      argument, e.g. ``shape=POLYGON((-10 50, 0 60, 10 50, 0 40, -10 50))``. See e.g.
+      https://en.wikipedia.org/wiki/Well-known_text for a description of the WKT format.
+
     .. note::
       Date/times are specified in the format: ``YYYY-MM-DDThh:mm:ss`` in which ``YYYY-MM-DD`` is a date and ``hh:mm:ss``
       is a time. A colon or space can be used instead of the 'T' separator (but if a space is used, the argument must be