Merge branch 'next' into release-0.08

Author: bekozi

doc/api.rst

@@ -310,6 +310,11 @@ Value             Description
 `False`           Maintain the :class:`~ocgis.RequestDataset`'s longitudinal domain.
 ================= ====================================================================================================
 
+add_auxiliary_files
+~~~~~~~~~~~~~~~~~~~
+
+If ``True``, create a new directory and add metadata and other informational files in addition to the converted file. If ``False``, write the target file only to :attr:`dir_output` and do not create a new directory.
+
 allow_empty
 ~~~~~~~~~~~
 

src/ocgis/api/operations.py

@@ -89,8 +89,8 @@ class OcgOperations(object):
     :param interpolate_spatial_bounds: If True and no bounds are available, attempt 
      to interpolate bounds from centroids.
     :type interpolate_spatial_bounds: bool
-    :param bool add_auxiliary_files: If True, create a new directory and add metadata 
-     and other informational files in addition to the converted file. If False, write
+    :param bool add_auxiliary_files: If ``True``, create a new directory and add metadata
+     and other informational files in addition to the converted file. If ``False``, write
      the target file only to :attr:`dir_output` and do not create a new directory.
     :param function callback: A function taking two parameters: ``percent_complete``
      and ``message``.

src/ocgis/api/subset.py

@@ -311,7 +311,7 @@ def _process_geometries_(self,itr,field,headers,value_keys,alias):
             ## if there is a slice, use it to subset the field.
             elif self.ops.slice is not None:
                 field = field.__getitem__(self.ops.slice)
-                
+
             ## see if the selection crs matches the field's crs
             if crs is not None and crs != field.spatial.crs:
                 geom = project_shapely_geometry(geom,crs.sr,field.spatial.crs.sr)
@@ -435,19 +435,22 @@ def _process_geometries_(self,itr,field,headers,value_keys,alias):
                     ## transform back to rotated pole if necessary
                     if original_rotated_pole_crs is not None:
                         if self.ops.output_crs is None and not isinstance(self.ops.output_crs,CFWGS84):
-                            ## we need to load the values before proceeding. source
-                            ## indices will disappear.
+                            # copy the spatial mask to the new spatial array
+                            spatial_mask_before_transform = deepcopy(sfield.spatial.get_mask())
+                            # need to load the values before proceeding. source indices will disappear.
                             for variable in sfield.variables.itervalues():
                                 variable.value
-                            ## reset the geometries
+                            # reset the geometries
                             sfield.spatial._geom = None
                             sfield.spatial.grid = get_rotated_pole_spatial_grid_dimension(
                              original_rotated_pole_crs,sfield.spatial.grid,inverse=True,
                              rc_original=original_row_column_metadata)
+                            # update the grid mask with the previous spatial mask
+                            sfield.spatial.grid.value.mask = spatial_mask_before_transform
                             ## update the uid mask to match the spatial mask
-                            sfield.spatial.uid = np.ma.array(sfield.spatial.uid,mask=sfield.spatial.get_mask())
+                            sfield.spatial.uid = np.ma.array(sfield.spatial.uid,mask=spatial_mask_before_transform)
                             sfield.spatial.crs = original_rotated_pole_crs
-                    
+
                     ## update the coordinate system of the data output
                     if self.ops.output_crs is not None:
                         ## if the geometry is not None, it may need to be projected to match

src/ocgis/conv/base.py

@@ -96,17 +96,22 @@ def write(self):
         ocgis_lh('starting write method',self._log,logging.DEBUG)
 
         unique_geometry_store = []
+
+        # indicates if user geometries should be written to file
+        write_ugeom = False
 
         try:
             build = True
-            if self._add_ugeom and self.ops is not None and self.ops.geom is not None:
-                write_geom = True
-            else:
-                write_geom = False
+
             for coll in iter(self.colls):
                 if build:
+
+                    # write the user geometries if configured and there is one present on the incoming collection.
+                    if self._add_ugeom and coll.geoms.values()[0] is not None:
+                        write_ugeom = True
+
                     f = self._build_(coll)
-                    if write_geom:
+                    if write_ugeom:
                         ugid_shp_name = self.prefix + '_ugid.shp'
                         ugid_csv_name = self.prefix + '_ugid.csv'
 
@@ -116,9 +121,14 @@ def write(self):
                         else:
                             fiona_path = os.path.join(self.outdir,ugid_shp_name)
                             csv_path = os.path.join(self.outdir,ugid_csv_name)
-                            
+
                         if coll.meta is None:
-                            fiona_properties = {'UGID':'int'}
+                            # convert the collection properties to fiona properties
+                            from fiona_ import FionaConverter
+                            fiona_properties = {}
+                            for k, v in coll.properties.values()[0].iteritems():
+                                fiona_properties[k] = FionaConverter.get_field_type(type(v))
+
                             fiona_schema = {'geometry':'MultiPolygon',
                                             'properties':fiona_properties}
                             fiona_meta = {'schema':fiona_schema,'driver':'ESRI Shapefile'}
@@ -148,7 +158,7 @@ def write(self):
 
                     build = False
                 self._write_coll_(f,coll)
-                if write_geom:
+                if write_ugeom:
                     ## write the overview geometries to disk
                     r_geom = coll.geoms.values()[0]
                     if isinstance(r_geom,Polygon):
@@ -187,7 +197,7 @@ def write(self):
                 ## this the exception we want to log
                 ocgis_lh(exc=e,logger=self._log)
             finally:
-                if write_geom:
+                if write_ugeom:
                     try:
                         fiona_object.close()
                     except UnboundLocalError:

src/ocgis/conv/fiona_.py

@@ -34,33 +34,44 @@ class FionaConverter(AbstractConverter):
                            np.int16:'int',
                            np.int32:'int',
                            str:'str'}
-    
+
+    @classmethod
+    def get_field_type(cls, the_type, key=None, fiona_conversion=None):
+        """
+        :param the_type: The target type object to map to a Fiona field type.
+        :type the_type: type object
+        :param key: The key to update the Fiona conversion map.
+        :type key: str
+        :param fiona_conversion: A dictionary used to convert Python values to Fiona-expected values.
+        :type fiona_conversion: dict
+        """
+
+        ret = None
+        for k, v in fiona.FIELD_TYPES_MAP.iteritems():
+            if the_type == v:
+                ret = k
+                break
+        if ret is None:
+            ret = cls._fiona_type_mapping[the_type]
+        if the_type in cls._fiona_conversion:
+            fiona_conversion.update({key.lower(): cls._fiona_conversion[the_type]})
+        return ret
+
     def _finalize_(self,f):
         f['fiona_object'].close()
 
     def _build_(self,coll):
         fiona_conversion = {}
 
-        def _get_field_type_(key,the_type):
-            ret = None
-            for k,v in fiona.FIELD_TYPES_MAP.iteritems():
-                if the_type == v:
-                    ret = k
-                    break
-            if ret is None:
-                ret = self._fiona_type_mapping[the_type]
-            if the_type in self._fiona_conversion:
-                fiona_conversion.update({key.lower():self._fiona_conversion[the_type]})
-            return(ret)
-        
         ## pull the fiona schema properties together by mapping fiona types to
         ## the data types of the first row of the output data file
         archetype_field = coll._archetype_field
         fiona_crs = archetype_field.spatial.crs.value
         geom,arch_row = coll.get_iter_dict().next()
         fiona_properties = OrderedDict()
         for header in coll.headers:
-            fiona_field_type = _get_field_type_(header,type(arch_row[header]))
+            fiona_field_type = self.get_field_type(type(arch_row[header]), key=header,
+                                                   fiona_conversion=fiona_conversion)
             fiona_properties.update({header.upper():fiona_field_type})
 
         ## we always want to convert the value. if the data is masked, it comes

src/ocgis/test/base.py

@@ -32,23 +32,28 @@ def _test_bin_dir(self):
         ret = os.path.join(base_dir, 'bin')
         return (ret)
 
-    def assertNumpyAll(self,arr1,arr2):
-        self.assertEqual(type(arr1),type(arr2))
-        if isinstance(arr1,np.ma.MaskedArray) or isinstance(arr2,np.ma.MaskedArray):
+    def assertNumpyAll(self, arr1, arr2):
+        self.assertEqual(type(arr1), type(arr2))
+        if isinstance(arr1, np.ma.MaskedArray) or isinstance(arr2, np.ma.MaskedArray):
             self.assertTrue(np.all(arr1.data == arr2.data))
             self.assertTrue(np.all(arr1.mask == arr2.mask))
-            self.assertEqual(arr1.fill_value,arr2.fill_value)
-            return(True)
+            self.assertEqual(arr1.fill_value, arr2.fill_value)
+            return True
         else:
-            return(self.assertTrue(np.all(arr1 == arr2)))
+            return self.assertTrue(np.all(arr1 == arr2))
 
     def assertNumpyAllClose(self, arr1, arr2):
         self.assertEqual(type(arr1), type(arr2))
         return self.assertTrue(np.allclose(arr1, arr2))
 
     def assertNumpyNotAll(self, arr1, arr2):
-        self.assertEqual(type(arr1), type(arr2))
-        return self.assertFalse(np.all(arr1 == arr2))
+        try:
+            self.assertNumpyAll(arr1, arr2)
+        except AssertionError:
+            ret = True
+        else:
+            raise AssertionError('Arrays are equivalent.')
+        return ret
 
     def assertDictEqual(self, d1, d2, msg=None):
         try:
@@ -163,7 +168,10 @@ def get_tdata():
         test_data.update(['misc', 'rotated_pole'], 'tas',
                          'tas_EUR-44_CCCma-CanESM2_rcp85_r1i1p1_SMHI-RCA4_v1_sem_209012-210011.nc',
                          key='rotated_pole_cccma')
-        return (test_data)
+        test_data.update(['misc', 'rotated_pole'], 'pr',
+                         'pr_EUR-11_CNRM-CERFACS-CNRM-CM5_historical_r1i1p1_CLMcom-CCLM4-8-17_v1_mon_198101-199012.nc',
+                         key='rotated_pole_cnrm_cerfacs')
+        return test_data
 
     def setUp(self):
         if self._reset_env: env.reset()

src/ocgis/test/test_ocgis/test_api/test_parms/test_definition.py

@@ -117,6 +117,33 @@ def test_geom_with_changing_select_ugid(self):
         g.select_ugid = [16,17]
         self.assertEqual(len(list(g.value)),2)
 
+    @staticmethod
+    def get_geometry_dictionaries():
+        coordinates = [('France', [2.8, 47.16]),
+                       ('Germany', [10.5, 51.29]),
+                       ('Italy', [12.2, 43.4])]
+        geom = []
+        for ugid, coordinate in enumerate(coordinates, start=1):
+            point = Point(coordinate[1][0], coordinate[1][1])
+            geom.append({'geom': point,
+                         'properties': {'UGID': ugid, 'COUNTRY': coordinate[0]}})
+        return geom
+
+    def test_geometry_dictionaries(self):
+        """Test geometry dictionaries come out appropriately once formatted."""
+
+        geom = self.get_geometry_dictionaries()
+        #todo: ensure geometry dictionaries have meta associated with them if more than a UGID is present in the properties
+        g = Geom(geom)
+
+        self.assertEqual(len(g.value), 3)
+
+        for gdict in g.value:
+            self.assertEqual(set(gdict.keys()), set(['crs', 'geom', 'properties']))
+            self.assertIsInstance(gdict['geom'], Point)
+            self.assertIsInstance(gdict['crs'], CFWGS84)
+            self.assertEqual(set(gdict['properties'].keys()), set(['UGID', 'COUNTRY']))
+
 
 class TestTimeRange(TestBase):
     _create_dir = False

src/ocgis/test/test_ocgis/test_api/test_subset.py

@@ -1,23 +1,40 @@
+from collections import OrderedDict
+from ocgis.conv.numpy_ import NumpyConverter
 from ocgis.test.base import TestBase
 import ocgis
 from ocgis.api.subset import SubsetOperation
 from ocgis.api.collection import SpatialCollection
 import itertools
+from ocgis.test.test_ocgis.test_api.test_parms.test_definition import TestGeom
 from ocgis.util.logging_ocgis import ProgressOcgOperations
 
 
-class Test(TestBase):
-    
+class TestSubsetOperation(TestBase):
+
     def get_operations(self):
         rd = self.test_data.get_rd('cancm4_tas')
-        slc = [None,[0,100],None,[0,10],[0,10]]
-        ops = ocgis.OcgOperations(dataset=rd,slice=slc)
-        return(ops)
-    
-    def test_constructor(self):
-        for rb,p in itertools.product([True,False],[None,ProgressOcgOperations()]):
-            sub = SubsetOperation(self.get_operations(),request_base_size_only=rb,
-                                  progress=p)
-            for ii,coll in enumerate(sub):
-                self.assertIsInstance(coll,SpatialCollection)
-        self.assertEqual(ii,0)
+        slc = [None, [0, 100], None, [0, 10], [0, 10]]
+        ops = ocgis.OcgOperations(dataset=rd, slice=slc)
+        return ops
+
+    def get_subset_operation(self):
+        geom = TestGeom.get_geometry_dictionaries()
+        rd = self.test_data.get_rd('cancm4_tas')
+        ops = ocgis.OcgOperations(dataset=rd, geom=geom, select_nearest=True)
+        subset = SubsetOperation(ops)
+        return subset
+
+    def test_init(self):
+        for rb, p in itertools.product([True, False], [None, ProgressOcgOperations()]):
+            sub = SubsetOperation(self.get_operations(), request_base_size_only=rb, progress=p)
+            for ii, coll in enumerate(sub):
+                self.assertIsInstance(coll, SpatialCollection)
+        self.assertEqual(ii, 0)
+
+    def test_geometry_dictionary(self):
+        """Test geometry dictionaries come out properly as collections."""
+
+        subset = self.get_subset_operation()
+        conv = NumpyConverter(subset, None, None)
+        coll = conv.write()
+        self.assertEqual(coll.properties, OrderedDict([(1, {'COUNTRY': 'France', 'UGID': 1}), (2, {'COUNTRY': 'Germany', 'UGID': 2}), (3, {'COUNTRY': 'Italy', 'UGID': 3})]))

src/ocgis/test/test_ocgis/test_conv/test_fiona_.py

@@ -0,0 +1,42 @@
+from collections import OrderedDict
+import os
+import fiona
+import ocgis
+from ocgis.api.subset import SubsetOperation
+from ocgis.conv.fiona_ import ShpConverter
+from ocgis.test.base import TestBase
+from ocgis.test.test_ocgis.test_api.test_parms.test_definition import TestGeom
+
+
+class TestShpConverter(TestBase):
+
+    def get_subset_operation(self):
+        geom = TestGeom.get_geometry_dictionaries()
+        rd = self.test_data.get_rd('cancm4_tas')
+        ops = ocgis.OcgOperations(dataset=rd, geom=geom, select_nearest=True, snippet=True)
+        subset = SubsetOperation(ops)
+        return subset
+
+    def test_attributes_copied(self):
+        """Test attributes in geometry dictionaries are properly accounted for in the converter."""
+
+        subset = self.get_subset_operation()
+        conv = ShpConverter(subset, self._test_dir, prefix='shpconv')
+        ret = conv.write()
+
+        path_ugid = os.path.join(self._test_dir, conv.prefix+'_ugid.shp')
+
+        with fiona.open(path_ugid) as source:
+            self.assertEqual(source.schema['properties'], OrderedDict([(u'COUNTRY', 'str'), (u'UGID', 'int:10')]))
+
+    def test_none_geom(self):
+        """Test a NoneType geometry will pass through the Fiona converter."""
+
+        rd = self.test_data.get_rd('cancm4_tas')
+        slc = [None, 0, None, [10, 20], [10, 20]]
+        ops = ocgis.OcgOperations(dataset=rd, slice=slc)
+        subset = SubsetOperation(ops)
+        conv = ShpConverter(subset, self._test_dir, prefix='shpconv')
+        ret = conv.write()
+        contents = os.listdir(self._test_dir)
+        self.assertEqual(len(contents), 5)