add function to select regions.

Author: gabrielaznar

climada/entity/exposures/base.py

@@ -16,7 +16,7 @@
 from climada.entity.tag import Tag
 from climada.util.coordinates import GridPoints
 from climada.util.config import CONFIG
-import climada.util.plot as plot
+import climada.util.plot as u_plot
 
 LOGGER = logging.getLogger(__name__)
 
@@ -26,7 +26,7 @@
            }
 """ Supported files format to read from """
 
-class Exposures(object):
+class Exposures():
     """Defines exposures attributes and basic methods. Loads from
     files with format defined in FILE_EXT.
 
@@ -149,10 +149,10 @@ def plot(self, ignore_zero=True, pop_name=True, buffer_deg=1.0,
             kwargs['reduce_C_function'] = np.sum
         if ignore_zero:
             pos_vals = self.value > 0
-            return plot.geo_bin_from_array(self.value[pos_vals], \
+            return u_plot.geo_bin_from_array(self.value[pos_vals], \
                 self.coord[pos_vals, :], cbar_label, title, pop_name, \
                 buffer_deg, extend, **kwargs)
-        return plot.geo_bin_from_array(self.value, self.coord, cbar_label, \
+        return u_plot.geo_bin_from_array(self.value, self.coord, cbar_label, \
             title, pop_name, buffer_deg, extend, **kwargs)
 
     def read(self, files, descriptions='', var_names=None):
@@ -261,6 +261,48 @@ def append(self, exposures):
             self.id[dup_id] = new_id
             new_id += 1
 
+    def select_region(self, reg_id):
+        """Return reference exposure of given region.
+
+        Parameters:
+            reg_id (int): region id to select
+
+        Returns:
+            Exposures
+        """
+        sel_idx = np.argwhere(self.region_id == reg_id)
+        if sel_idx.size:
+            sel_idx = sel_idx.reshape((sel_idx.size,))
+        else:
+            LOGGER.info('No exposure with region id %s.', reg_id)
+            return None
+
+        sel_exp = Exposures()
+        sel_exp.tag = self.tag
+        sel_exp.description = 'Region ' + str(reg_id)
+        sel_exp.ref_year = self.ref_year
+        sel_exp.value_unit = self.value_unit
+        # Obligatory variables
+        sel_exp.coord = self.coord[sel_idx, :]
+        sel_exp.value = self.value[sel_idx]
+        sel_exp.impact_id = self.impact_id[sel_idx]
+        sel_exp.id = self.id[sel_idx]
+        # Optional variables.
+        if self.deductible.size:
+            sel_exp.deductible = self.deductible[sel_idx]
+        if self.cover.size:
+            sel_exp.cover = self.cover[sel_idx]
+        if self.category_id.size:
+            sel_exp.category_id = self.category_id[sel_idx]
+        if self.region_id.size:
+            sel_exp.region_id = self.region_id[sel_idx]
+
+        sel_exp.assigned = dict()
+        for key, value in self.assigned.items():
+            sel_exp.assigned[key] = value[sel_idx]
+
+        return sel_exp
+
     @staticmethod
     def get_sup_file_format():
         """ Get supported file extensions that can be read.

climada/entity/exposures/black_marble.py

@@ -15,6 +15,7 @@
 import requests
 import shapely.vectorized
 from cartopy.io import shapereader
+from iso3166 import countries as iso_cntry
 
 from climada.entity.exposures.base import Exposures
 from climada.util.files_handler import download_file
@@ -59,15 +60,19 @@
 
 class BlackMarble(Exposures):
     """Defines exposures from night light intensity, GDP and income group.
+    Attribute region_id is defined as:
+    - United Nations Statistics Division (UNSD) 3-digit equivalent numeric code
+    - 0 if country not found in UNSD.
+    - -1 for water
     """
 
     def __init__(self):
         """ Empty initializer. """
         Exposures.__init__(self)
 
-    def set_countries(self, countries, \
-        ref_year=CONFIG['entity']['present_ref_year'], res_km=1, \
-        sea_res=(0, 1), from_hr=None, **kwargs):
+    def set_countries(self, countries,
+                      ref_year=CONFIG['entity']['present_ref_year'],
+                      res_km=1, sea_res=(0, 1), from_hr=None, **kwargs):
         """ Model countries using values at reference year. If GDP or income
         group not available for that year, consider the value of the closest
         available year.
@@ -80,10 +85,11 @@ def set_countries(self, countries, \
             res_km (float, optional): approx resolution in km. Default: 1km.
             sea_res (tuple, optional): (sea_coast_km, sea_res_km), where first
                 parameter is distance from coast to fill with water and second
-                parameter is resolution between sea points
-            from_hr (bool, optional)
+                parameter is resolution between sea points.
+            from_hr (bool, optional): force to use higher resolution image,
+                independently of its year of acquisition.
             kwargs (optional): 'gdp' and 'inc_grp' dictionaries with keys the
-                country ISO_alpha3 code. If provided, these are used
+                country ISO_alpha3 code. If provided, these are used.
         """
         self.clear()
         shp_file = shapereader.natural_earth(resolution='10m',
@@ -97,23 +103,38 @@ def set_countries(self, countries, \
         nightlight, coord_nl, fn_nl, res_fact = get_nightlight(ref_year, \
             cntry_info, res_km, from_hr)
 
-        # TODO parallize per thread?
         for cntry_iso, cntry_val in cntry_info.items():
             LOGGER.info('Processing country %s.', cntry_val[1])
             self.append(self._set_one_country(cntry_val, nightlight, \
                 coord_nl, fn_nl, res_fact, res_km, cntry_admin1[cntry_iso]))
 
         add_sea(self, sea_res)
 
-    def set_region(self, centroids, ref_year=2013, resolution=1):
+    def set_region(self, centroids,
+                   ref_year=CONFIG['entity']['present_ref_year'], res_km=1):
         """ Model a specific region given by centroids."""
         # TODO: accept input centroids as well
         raise NotImplementedError
 
-    def fill_centroids(self):
-        """ From coordinates information, generate Centroids instance."""
-        # add sea in lower resolution
-        raise NotImplementedError
+    def select_region(self, reg_id):
+        """ Select exposures with input region.
+
+        Parameters:
+            reg_id (int, str): integer iso equivalent country numeric code or
+                string iso alpha-3 or alpha-2 code or country name.
+
+        Returns:
+            Exposures
+        """
+        if isinstance(reg_id, int):
+            return Exposures.select_region(self, reg_id)
+
+        try:
+            return Exposures.select_region(self, \
+                int(iso_cntry.get(reg_id).numeric))
+        except KeyError:
+            LOGGER.info('No country %s found.', reg_id)
+            return None
 
     @staticmethod
     def _set_one_country(cntry_info, nightlight, coord_nl, fn_nl, res_fact,
@@ -155,8 +176,9 @@ def _set_one_country(cntry_info, nightlight, coord_nl, fn_nl, res_fact,
         return exp_bkmrb
 
 def country_iso_geom(countries, shp_file):
-    """ Get country ISO alpha_3, country id (defined as country appearance
-    order in natural earth shape file) and country's geometry.
+    """ Get country ISO alpha_3, country id (defined as the United Nations
+    Statistics Division (UNSD) 3-digit equivalent numeric codes and 0 if
+    country not found) and country's geometry shape.
 
     Parameters:
         countries (list or dict): list of country names (admin0) or dict
@@ -172,8 +194,6 @@ def country_iso_geom(countries, shp_file):
     for info_idx, info in enumerate(list_records):
         countries_shp[info.attributes['ADMIN'].title()] = info_idx
 
-    num_codes = [iso3 for iso3 in wb.country_codes if len(iso3) == 3]
-
     cntry_info = dict()
     cntry_admin1 = dict()
     if isinstance(countries, list):
@@ -198,10 +218,10 @@ def country_iso_geom(countries, shp_file):
             raise ValueError
         iso3 = list_records[country_idx].attributes['ADM0_A3']
         try:
-            iso_num = num_codes.index(iso3)
-        except ValueError:
-            iso_num = len(num_codes)
-        cntry_info[iso3] = [iso_num, country_name.title(),
+            cntry_id = int(iso_cntry.get(iso3).numeric)
+        except KeyError:
+            cntry_id = 0
+        cntry_info[iso3] = [cntry_id, country_name.title(),
                             list_records[country_idx].geometry]
         cntry_admin1[iso3] = _fill_admin1_geom(iso3, admin1_rec, prov_list)
 
@@ -261,7 +281,7 @@ def get_nightlight(ref_year, cntry_info, res_km, from_hr=None):
         else:
             nl_year = 2012
         LOGGER.info("Nightlights from NASA's earth observatory for " + \
-                    "year %s.", nl_year)
+                    "year %s.", str(nl_year))
         res_fact = DEF_RES_NASA_KM/res_km
         geom = [info[2] for info in cntry_info.values()]
         geom = shapely.ops.cascaded_union(geom)
@@ -281,7 +301,7 @@ def get_nightlight(ref_year, cntry_info, res_km, from_hr=None):
         elif ref_year > 2013:
             nl_year = 2013
         LOGGER.info("Nightlights from NOAA's earth observation group "+ \
-                    "for year %s.", nl_year)
+                    "for year %s.", str(nl_year))
         res_fact = DEF_RES_NOAA_KM/res_km
         # nightlight intensity with 30 arcsec resolution
         nightlight, coord_nl, fn_nl = nl_utils.load_nightlight_noaa(nl_year)
@@ -301,7 +321,7 @@ def add_sea(exp, sea_res):
         return
 
     LOGGER.info("Adding sea at %s km resolution and %s km distance from " + \
-        "coast.", sea_res[1], sea_res[0])
+        "coast.", str(sea_res[1]), str(sea_res[0]))
 
     sea_res = (sea_res[0]/ONE_LAT_KM, sea_res[1]/ONE_LAT_KM)
 

climada/entity/exposures/test/test_base.py

@@ -347,11 +347,44 @@ def test_check_wrongId_fail(self):
         self.assertIn('There are exposures with the same identifier.',
                       cm.output[0])
 
+class TestSelectRegion(unittest.TestCase):
+    """Test select_region from the Exposures class"""
+    def test_sel_reg_pass(self):
+        """Select region"""
+        expo = good_exposures()
+        sel_expo = expo.select_region(1)
+
+        self.assertEqual(sel_expo.value.size, 1)
+        self.assertEqual(sel_expo.value[0], 1)
+        self.assertEqual(sel_expo.id.size, 1)
+        self.assertEqual(sel_expo.id[0], 1)
+        self.assertEqual(sel_expo.impact_id.size, 1)
+        self.assertEqual(sel_expo.impact_id[0], 1)
+        self.assertEqual(sel_expo.deductible.size, 1)
+        self.assertEqual(sel_expo.deductible[0], 1)
+        self.assertEqual(sel_expo.category_id.size, 1)
+        self.assertEqual(sel_expo.category_id[0], 1)
+        self.assertEqual(sel_expo.region_id.size, 1)
+        self.assertEqual(sel_expo.region_id[0], 1)
+        self.assertEqual(sel_expo.cover.size, 0)
+        self.assertEqual(sel_expo.assigned['TC'].size, 1)
+        self.assertEqual(sel_expo.assigned['TC'][0], 1)
+        self.assertEqual(sel_expo.coord.shape[0], 1)
+        self.assertEqual(sel_expo.coord[0, 0], 1)
+        self.assertEqual(sel_expo.coord[0, 1], 2)
+
+    def test_sel_wrong_pass(self):
+        """Select non-existent region"""
+        expo = good_exposures()
+        sel_expo = expo.select_region(5)
+        self.assertEqual(sel_expo, None)
+
 # Execute Tests
 TESTS = unittest.TestLoader().loadTestsFromTestCase(TestChecker)
 TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAssign))
 TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestAppend))
 TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestReadParallel))
 TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestConstructor))
 TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestRemove))
+TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestSelectRegion))
 unittest.TextTestRunner(verbosity=2).run(TESTS)

climada/entity/exposures/test/test_black_marble.py

@@ -35,10 +35,10 @@ def test_che_kos_pass(self):
         self.assertEqual(len(iso_name), len(country_name))
         self.assertTrue('CHE'in iso_name)
         self.assertTrue('KOS'in iso_name)
-        self.assertEqual(iso_name['CHE'][0], 41)
+        self.assertEqual(iso_name['CHE'][0], 756)
         self.assertEqual(iso_name['CHE'][1], 'Switzerland')
         self.assertIsInstance(iso_name['CHE'][2], shapely.geometry.multipolygon.MultiPolygon)
-        self.assertEqual(iso_name['KOS'][0], 252)
+        self.assertEqual(iso_name['KOS'][0], 0)
         self.assertEqual(iso_name['KOS'][1], 'Kosovo')
         self.assertIsInstance(iso_name['KOS'][2], shapely.geometry.multipolygon.MultiPolygon)
 
@@ -48,7 +48,7 @@ def test_haiti_pass(self):
         iso_name, _ = country_iso_geom(country_name, SHP_FILE)
 
         self.assertEqual(len(iso_name), len(country_name))
-        self.assertEqual(iso_name['HTI'][0], 100)
+        self.assertEqual(iso_name['HTI'][0], 332)
         self.assertEqual(iso_name['HTI'][1], 'Haiti')
         self.assertIsInstance(iso_name['HTI'][2], shapely.geometry.multipolygon.MultiPolygon)
 
@@ -64,7 +64,7 @@ def test_bolivia_pass(self):
         iso_name, _ = country_iso_geom(country_name, SHP_FILE)
 
         self.assertEqual(len(iso_name), len(country_name))
-        self.assertEqual(iso_name['BOL'][0], 31)
+        self.assertEqual(iso_name['BOL'][0], 68)
         self.assertEqual(iso_name['BOL'][1], 'Bolivia')
         self.assertIsInstance(iso_name['BOL'][2], shapely.geometry.multipolygon.MultiPolygon)
 
@@ -73,6 +73,62 @@ def test_korea_pass(self):
         country_name = ['Korea']
         with self.assertRaises(ValueError):
             country_iso_geom(country_name, SHP_FILE)
+    
+    def test_select_reg_pass(self):
+        """Select counrty """
+        exp = BlackMarble()
+        exp.value = np.arange(0, 20, 0.1)
+        exp.coord = np.ones((200, 2))
+        exp.impact_id = np.ones((200,))
+        exp.id = np.arange(200)
+        
+        exp.region_id = np.ones((200,))
+        exp.region_id[:101] = exp.region_id[:101]*780
+        exp.region_id[101:] = exp.region_id[101:]*894
+        
+        sel_exp = exp.select_region('TTO')
+        self.assertEqual(sel_exp.value.size, 101)
+        sel_exp = exp.select_region('Trinidad and Tobago')
+        self.assertEqual(sel_exp.value.size, 101)
+        sel_exp = exp.select_region('TT')
+        self.assertEqual(sel_exp.value.size, 101)
+
+        sel_exp = exp.select_region('ZMB')
+        self.assertEqual(sel_exp.value.size, 99)
+        sel_exp = exp.select_region('Zambia')
+        self.assertEqual(sel_exp.value.size, 99)
+        sel_exp = exp.select_region('ZM')
+        self.assertEqual(sel_exp.value.size, 99)
+
+    def test_select_not_included_pass(self):
+        """Select country that is not contained."""
+        exp = BlackMarble()
+        exp.value = np.arange(0, 20, 0.1)
+        exp.coord = np.ones((200, 2))
+        exp.impact_id = np.ones((200,))
+        exp.id = np.arange(200)
+        
+        exp.region_id = np.ones((200,))
+        exp.region_id[:101] = exp.region_id[:101]*780
+        exp.region_id[101:] = exp.region_id[101:]*894
+        
+        sel_exp = exp.select_region('UGA')
+        self.assertEqual(sel_exp, None)
+
+    def test_select_wrong_fail(self):
+        """Select country that is not contained."""
+        exp = BlackMarble()
+        exp.value = np.arange(0, 20, 0.1)
+        exp.coord = np.ones((200, 2))
+        exp.impact_id = np.ones((200,))
+        exp.id = np.arange(200)
+        
+        exp.region_id = np.ones((200,))
+        exp.region_id[:101] = exp.region_id[:101]*780
+        exp.region_id[101:] = exp.region_id[101:]*894
+        
+        sel_exp = exp.select_region('TTU')
+        self.assertEqual(sel_exp, None)
 
 class TestProvinces(unittest.TestCase):
     """Tst black marble with admin1."""
@@ -371,8 +427,9 @@ def test_set_econ_indicators_pass(self):
         self.assertAlmostEqual(exp.value.sum(), gdp*(inc_grp+1), 5)
 
 # Execute Tests
-TESTS = unittest.TestLoader().loadTestsFromTestCase(TestEconIndices)
-TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestCountryIso))
-TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestNightLight))
-TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestProvinces))
+TESTS = unittest.TestLoader().loadTestsFromTestCase(TestCountryIso)
+#TESTS = unittest.TestLoader().loadTestsFromTestCase(TestEconIndices)
+#TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestCountryIso))
+#TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestNightLight))
+#TESTS.addTests(unittest.TestLoader().loadTestsFromTestCase(TestProvinces))
 unittest.TextTestRunner(verbosity=2).run(TESTS)