define meta attribute in exposures

gabrielaznar · gabrielaznar · commit 7804bd9ffc63 · 2019-06-18T20:18:55.000+02:00
diff --git a/climada/entity/exposures/base.py b/climada/entity/exposures/base.py
@@ -87,6 +87,10 @@ class Exposures(GeoDataFrame):
             If not provided, set to default 'if_' with ids 1 in check().
         geometry (pd.Series, optional): geometry of type Point of each instance.
             Computed in method set_geometry_points().
+        meta (dict): dictionary containing corresponding raster properties (if any):
+            width, height, crs and transform must be present at least (transform needs
+            to contain upper left corner!). Exposures might not contain all the points
+            of the corresponding raster.
         deductible (pd.Series, optional): deductible value for each exposure
         cover (pd.Series, optional): cover value for each exposure
         category_id (pd.Series, optional): category id for each exposure
@@ -95,7 +99,8 @@ class Exposures(GeoDataFrame):
             TC. There might be different hazards defined: centr_TC, centr_FL, ...
             Computed in method assign_centroids().
     """
-    _metadata = GeoDataFrame._metadata + ['tag', 'ref_year', 'value_unit']
+    _metadata = GeoDataFrame._metadata + ['tag', 'ref_year', 'value_unit',
+                                          'meta']
 
     vars_oblig = ['value', 'latitude', 'longitude']
     """Name of the variables needed to compute the impact."""
@@ -139,6 +144,8 @@ def check(self):
                     self.ref_year = DEF_REF_YEAR
                 elif var == 'value_unit':
                     self.value_unit = DEF_VALUE_UNIT
+                elif var == 'meta':
+                    self.meta = None
                 LOGGER.info('%s metadata set to default value: %s', var, self.__dict__[var])
 
         for var in self.vars_oblig:
@@ -215,14 +222,14 @@ def assign_centroids(self, hazard, method='NN', distance='haversine',
         self[INDICATOR_CENTR + hazard.tag.haz_type] = assigned
 
     def set_geometry_points(self, scheduler=None):
-        """ Set geometry attribute of GeoDataFrame from latitude and longitude
-        attributes.
+        """ Set geometry attribute of GeoDataFrame with Points from latitude and
+        longitude attributes.
 
         Parameter:
             scheduler (str): used for dask map_partitions. “threads”,
                 “synchronous” or “processes”
         """
-        LOGGER.info('Setting geometry attribute.')
+        LOGGER.info('Setting geometry points.')
         def apply_point(df_exp):
             return df_exp.apply((lambda row: Point(row.longitude, row.latitude)), axis=1)
         if not scheduler:
@@ -244,7 +251,7 @@ def set_lat_lon(self):
     def set_from_raster(self, file_name, band=1, src_crs=None, window=False,
                         geometry=False, dst_crs=False, transform=None,
                         width=None, height=None, resampling=Resampling.nearest):
-        """ Read raster data and set latitude, longitude and value
+        """ Read raster data and set latitude, longitude, value and meta
 
         Parameters:
             file_name (str): file name containing values
@@ -259,9 +266,6 @@ def set_from_raster(self, file_name, band=1, src_crs=None, window=False,
             height (float): number of lats for transform
             resampling (rasterio.warp,.Resampling optional): resampling
                 function used for reprojection to dst_crs
-
-        Returns:
-            dict (meta raster information)
         """
         self.__init__()
         self.tag = Tag()
@@ -281,7 +285,7 @@ def set_from_raster(self, file_name, band=1, src_crs=None, window=False,
         self['longitude'] = x_grid.flatten()
         self['latitude'] = y_grid.flatten()
         self['value'] = value.reshape(-1)
-        return meta
+        self.meta = meta
 
     def plot_scatter(self, mask=None, ignore_zero=False, pop_name=True,
                      buffer=0.0, extend='neither', **kwargs):
diff --git a/climada/entity/exposures/black_marble.py b/climada/entity/exposures/black_marble.py
@@ -34,9 +34,10 @@
 
 from climada.entity.tag import Tag
 from climada.entity.exposures.base import Exposures, INDICATOR_IF
-from climada.util.constants import SYSTEM_DIR
+from climada.util.constants import SYSTEM_DIR, DEF_CRS
 from climada.entity.exposures import nightlight as nl_utils
 from climada.util.finance import gdp, income_group
+from climada.util.coordinates import points_to_raster
 
 LOGGER = logging.getLogger(__name__)
 
@@ -104,15 +105,18 @@ def set_countries(self, countries, ref_year=2016, res_km=None, from_hr=None,
             tag.description += ("{} {:d} GDP: {:.3e} income group: {:d} \n").\
                 format(cntry_val[1], cntry_val[3], cntry_val[4], cntry_val[5])
 
-        Exposures.__init__(self, gpd.GeoDataFrame(pd.concat(bkmrbl_list,
-                                                            ignore_index=True)))
+        Exposures.__init__(self, gpd.GeoDataFrame(pd.concat(bkmrbl_list, \
+            ignore_index=True)), crs=DEF_CRS)
 
         # set metadata
         self.ref_year = ref_year
         self.tag = tag
         self.tag.file_name = fn_nl
         self.value_unit = 'USD'
-        self.crs = {'init': 'epsg:4326'}
+        rows, cols, ras_trans = points_to_raster((self.longitude.min(), \
+            self.latitude.min(), self.longitude.max(), self.latitude.max()), \
+            coord_nl[0, 1])
+        self.meta = {'width':cols, 'height':rows, 'crs':self.crs, 'transform':ras_trans}
 
     @staticmethod
     def _set_one_country(cntry_info, nightlight, coord_nl, res_fact,
diff --git a/climada/entity/exposures/litpop.py b/climada/entity/exposures/litpop.py
@@ -38,7 +38,8 @@
 from climada.entity.exposures.base import Exposures, INDICATOR_IF
 from climada.entity.exposures import gpw_import
 from climada.util.finance import gdp, income_group, wealth2gdp, world_bank_wealth_account
-from climada.util.constants import SYSTEM_DIR
+from climada.util.constants import SYSTEM_DIR, DEF_CRS
+from climada.util.coordinates import points_to_raster, get_resolution
 
 logging.root.setLevel(logging.DEBUG)
 LOGGER = logging.getLogger(__name__)
@@ -261,11 +262,15 @@ def set_country(self, countries, **args):
                 min(GPW_YEARS, key=lambda x:abs(x-reference_year)), \
                 min(BM_YEARS, key=lambda x:abs(x-reference_year)), \
                 exponents[0], exponents[1])
-        Exposures.__init__(self, gpd.GeoDataFrame(pd.concat(lp_cntry,
-                                                            ignore_index=True)))
+        Exposures.__init__(self, gpd.GeoDataFrame(pd.concat(lp_cntry, \
+            ignore_index=True)), crs=DEF_CRS)
         self.ref_year = reference_year
         self.tag = tag
         self.value_unit = 'USD'
+        rows, cols, ras_trans = points_to_raster((self.longitude.min(), \
+            self.latitude.min(), self.longitude.max(), self.latitude.max()), \
+            min(get_resolution(self.latitude, self.longitude)))
+        self.meta = {'width':cols, 'height':rows, 'crs':self.crs, 'transform':ras_trans}
 
         if check_plot == 1:
             self.plot_log(admin1_plot=0)
@@ -1875,11 +1880,11 @@ def exposure_set_admin1(exposure):
     """
     Input:
         exposure: exposure instance
-        
+
     Returns:
         exposure: exposure instance with 2 extra columns: admin1 & admin1_ID
     """
-    
+
     exposure['admin1'] = pd.Series()
     exposure['admin1_ID'] = pd.Series()
     count = 0
diff --git a/climada/test/test_LitPop.py b/climada/test/test_LitPop.py
@@ -25,6 +25,7 @@
 from climada.entity.exposures import litpop as lp
 from climada.entity.exposures import gpw_import
 from climada.util.finance import world_bank_wealth_account
+from climada.util.coordinates import equal_crs
 
 # ---------------------
 class TestLitPopExposure(unittest.TestCase):
@@ -48,6 +49,16 @@ def test_switzerland300_pass(self):
         self.assertIn('GPW-year=2015', ent.tag.description)
         self.assertIn('BM-year=2016', ent.tag.description)
         self.assertIn('exp=[1, 1]', ent.tag.description)
+        self.assertTrue(equal_crs(ent.crs['init'], {'init': 'epsg:4326'}))
+        self.assertEqual(ent.meta['width'], 54)
+        self.assertEqual(ent.meta['height'], 23)
+        self.assertTrue(equal_crs(ent.meta['crs'], {'init': 'epsg:4326'}))
+        self.assertAlmostEqual(ent.meta['transform'][0], 0.08333333333333333)
+        self.assertAlmostEqual(ent.meta['transform'][1], 0)
+        self.assertAlmostEqual(ent.meta['transform'][2], 5.9166666666666)
+        self.assertAlmostEqual(ent.meta['transform'][3], 0)
+        self.assertAlmostEqual(ent.meta['transform'][4], -0.08333333333333333)
+        self.assertAlmostEqual(ent.meta['transform'][5], 47.7499999999999)
 
     def test_switzerland30normPop_pass(self):
         """Create LitPop entity for Switzerland on 30 arcsec:"""
diff --git a/climada/test/test_blackmarble.py b/climada/test/test_blackmarble.py
@@ -28,6 +28,7 @@
 from climada.entity.exposures.nightlight import load_nightlight_nasa, \
 load_nightlight_noaa, NOAA_BORDER, cut_nl_nasa
 from climada.entity.exposures import nightlight as nl_utils
+from climada.util.coordinates import equal_crs
 
 class Test2013(unittest.TestCase):
     """Test black marble of previous in 2013."""
@@ -46,7 +47,16 @@ def test_spain_pass(self):
         self.assertIn("Processing country Spain.", cm.output[1])
         self.assertIn("Generating resolution of approx 1 km.", cm.output[2])
         self.assertTrue(np.isclose(ent.value.sum(), 1.362e+12*(4+1), 4))
-        self.assertEqual(ent.crs, {'init': 'epsg:4326'})
+        self.assertTrue(equal_crs(ent.crs['init'], {'init': 'epsg:4326'}))
+        self.assertEqual(ent.meta['width'], 2699)
+        self.assertEqual(ent.meta['height'], 1938)
+        self.assertTrue(equal_crs(ent.meta['crs'], {'init': 'epsg:4326'}))
+        self.assertAlmostEqual(ent.meta['transform'][0], 0.008333333333333333)
+        self.assertAlmostEqual(ent.meta['transform'][1], 0)
+        self.assertAlmostEqual(ent.meta['transform'][2], -18.1625000000000)
+        self.assertAlmostEqual(ent.meta['transform'][3], 0)
+        self.assertAlmostEqual(ent.meta['transform'][4], -0.008333333333333333)
+        self.assertAlmostEqual(ent.meta['transform'][5], 43.79583333333333)
 
     def test_sint_maarten_pass(self):
         country_name = ['Sint Maarten']
@@ -55,15 +65,15 @@ def test_sint_maarten_pass(self):
             ent.set_countries(country_name, 2013, res_km=0.2)
         self.assertIn('GDP SXM 2014: 3.658e+08.', cm.output[0])
         self.assertIn('Income group SXM 2013: 4.', cm.output[1])
-        self.assertEqual(ent.crs, {'init': 'epsg:4326'})
+        self.assertTrue(equal_crs(ent.crs['init'], {'init': 'epsg:4326'}))
 
         with self.assertLogs('climada.entity.exposures.black_marble', level='INFO') as cm:
             ent.set_countries(country_name, 2013, res_km=0.2)
         self.assertIn("Nightlights from NOAA's earth observation group for year 2013.", cm.output[0])
         self.assertIn("Processing country Sint Maarten.", cm.output[1])
         self.assertIn("Generating resolution of approx 0.2 km.", cm.output[2])
         self.assertAlmostEqual(ent.value.sum(), 3.658e+08*(4+1))
-        self.assertEqual(ent.crs, {'init': 'epsg:4326'})
+        self.assertTrue(equal_crs(ent.crs['init'], {'init': 'epsg:4326'}))
 
     def test_anguilla_pass(self):
         country_name = ['Anguilla']
@@ -72,7 +82,7 @@ def test_anguilla_pass(self):
         self.assertEqual(ent.ref_year, 2013)
         self.assertIn("Anguilla 2013 GDP: 1.754e+08 income group: 3", ent.tag.description)
         self.assertAlmostEqual(ent.value.sum(), 1.754e+08*(3+1))
-        self.assertEqual(ent.crs, {'init': 'epsg:4326'})
+        self.assertTrue(equal_crs(ent.crs['init'], {'init': 'epsg:4326'}))
 
 class Test1968(unittest.TestCase):
     """Test black marble of previous years to 1992."""
@@ -90,7 +100,7 @@ def test_switzerland_pass(self):
         self.assertTrue("Processing country Switzerland." in cm.output[-2])
         self.assertTrue("Generating resolution of approx 0.5 km." in cm.output[-1])
         self.assertTrue(np.isclose(ent.value.sum(), 1.894e+10*(4+1), 4))
-        self.assertEqual(ent.crs, {'init': 'epsg:4326'})
+        self.assertTrue(equal_crs(ent.crs['init'], {'init': 'epsg:4326'}))
 
 class Test2012(unittest.TestCase):
     """Test year 2012 flags."""
@@ -126,7 +136,7 @@ def test_from_hr_flag_pass(self):
         self.assertTrue(np.isclose(ent.value.sum(), 8.740e+11*(3+1), 4))
         size3 = ent.value.size
         self.assertEqual(size1, size3)
-        self.assertEqual(ent.crs, {'init': 'epsg:4326'})
+        self.assertTrue(equal_crs(ent.crs['init'], {'init': 'epsg:4326'}))
 
 class BMFuncs(unittest.TestCase):
     """Test plot functions."""
