Auto assign basin projection attributes and basin center

requirements.txt

@@ -7,7 +7,7 @@ numpy>=1.14.0
 pandas>=0.22.0
 pytz==2017.3
 scipy>=1.0.0
-utm==0.4.0
+utm==0.5.0
 git+https://github.com/USDA-ARS-NWRC/PyKrige.git@master
 spatialnc>=0.2.12,<0.3.0
 weather-forecast-retrieval

smrf/data/loadGrid.py

@@ -51,10 +51,6 @@ def __init__(self, dataConfig, topo, start_date, end_date,
         # degree offset for a buffer around the model domain
         self.offset = 0.1
 
-        self.force_zone_number = None
-        if 'zone_number' in dataConfig:
-            self.force_zone_number = dataConfig['zone_number']
-
         # The data that will be output
         self.variables = ['air_temp', 'vapor_pressure', 'precip', 'wind_speed',
                           'wind_direction', 'cloud_factor', 'thermal']
@@ -63,18 +59,16 @@ def __init__(self, dataConfig, topo, start_date, end_date,
         # the model domain are used
         self.x = topo.x
         self.y = topo.y
-        self.lat = topo.topoConfig['basin_lat']
-        self.lon = topo.topoConfig['basin_lon']
 
-        # get the zone number and the bounding box
-        u = utm.from_latlon(topo.topoConfig['basin_lat'],
-                            topo.topoConfig['basin_lon'],
-                            self.force_zone_number)
-        self.zone_number = u[2]
-        self.zone_letter = u[3]
+        self.lat = topo.basin_lat
+        self.lon = topo.basin_long
+        self.zone_number = topo.zone_number
+        self.northern = topo.northern_hemisphere
 
-        ur = np.array(utm.to_latlon(np.max(self.x), np.max(self.y), self.zone_number, self.zone_letter))
-        ll = np.array(utm.to_latlon(np.min(self.x), np.min(self.y), self.zone_number, self.zone_letter))
+        ur = np.array(utm.to_latlon(np.max(self.x), np.max(self.y),
+                                      self.zone_number, northern=self.northern))
+        ll = np.array(utm.to_latlon(np.min(self.x), np.min(self.y),
+                                      self.zone_number, northern=self.northern))
 
         buff = 0.1 # buffer of bounding box in degrees
         ur += buff
@@ -98,11 +92,12 @@ def model_domain_grid(self):
         dlat = np.zeros((2,))
         dlon = np.zeros_like(dlat)
         dlat[0], dlon[0] = utm.to_latlon(np.min(self.x), np.min(self.y),
-                                         int(self.dataConfig['zone_number']),
-                                         self.dataConfig['zone_letter'])
+                                         self.zone_number,
+                                         northern=self.northern)
+
         dlat[1], dlon[1] = utm.to_latlon(np.max(self.x), np.max(self.y),
-                                         int(self.dataConfig['zone_number']),
-                                         self.dataConfig['zone_letter'])
+                                         self.zone_number,
+                                         northern=self.northern)
         # add a buffer
         dlat[0] -= self.offset
         dlat[1] += self.offset
@@ -145,7 +140,7 @@ def load_from_hrrr(self):
             self.end_date,
             self.bbox,
             output_dir=self.dataConfig['hrrr_directory'],
-            force_zone_number=self.force_zone_number,
+            force_zone_number=self.zone_number,
             forecast=fcast,
             forecast_flag=self.forecast_flag,
             day_hour=self.day_hour)
@@ -249,7 +244,7 @@ def load_from_netcdf(self):
         metadata['longitude'] = mlon.flatten()
         metadata['elevation'] = mhgt.flatten()
         metadata = metadata.apply(apply_utm,
-                                  args=(self.force_zone_number,),
+                                  args=(self.zone_number,),
                                   axis=1)
 
         self.metadata = metadata
@@ -358,7 +353,7 @@ def load_from_wrf(self):
         metadata['longitude'] = mlon.flatten()
         metadata['elevation'] = mhgt.flatten()
         metadata = metadata.apply(apply_utm,
-                                  args=(self.force_zone_number,),
+                                  args=(self.zone_number,),
                                   axis=1)
 
         self.metadata = metadata

smrf/data/loadTopo.py

@@ -7,9 +7,10 @@
 import numpy as np
 from netCDF4 import Dataset
 from spatialnc import ipw
+from utm import to_latlon
 
 from smrf.utils import gradient
-
+from smrf.utils.utils import get_center
 
 class topo():
     """
@@ -109,6 +110,27 @@ def readNetCDF(self):
         self.y = f.variables['y'][:]
         [self.X, self.Y] = np.meshgrid(self.x, self.y)
 
+
+        # Calculate the center of the basin
+        self.cx, self.cy  = get_center(f, mask_name='mask')
+
+        # Is the modeling domain in the northern hemisphere
+        self.northern_hemisphere = self.topoConfig['northern_hemisphere']
+
+        # Assign the UTM zone
+        self.zone_number = int(f.variables['projection'].utm_zone_number)
+
+        # Calculate the lat long
+        self.basin_lat, self.basin_long = to_latlon(self.cx,
+                                            self.cy,
+                                            self.zone_number,
+                                            northern=self.northern_hemisphere)
+
+        self._logger.info('Domain center in UTM Zone {:d} = {:0.1f}m, {:0.1f}m'
+                          ''.format(self.zone_number, self.cx, self.cy))
+        self._logger.info('Domain center as Latitude/Longitude = {:0.5f}, '
+                          '{:0.5f}'.format(self.basin_lat, self.basin_long))
+
         f.close()
 
     def stoporadInput(self):
@@ -185,7 +207,7 @@ def _gradient(self, demFile, gradientFile):
             raise OSError('gradient failed')
 
     def gradient(self, gfile):
-        """ 
+        """
         Calculate the gradient and aspect
 
         Args:

smrf/framework/CoreConfig.ini

@@ -19,16 +19,6 @@
 
 [topo]
 
-# TODO remove and auto calculate #97
-basin_lat:
-type = float,
-description = Latitude of the center of the basin used for sun angle calcs.
-
-# TODO remove and auto calculate #97
-basin_lon:
-type = float,
-description = Longitude of the center of the basin used for sun angle calcs.
-
 filename:
 type = CriticalFilename,
 description = A netCDF file containing all veg info and dem.
@@ -37,10 +27,15 @@ gradient_method:
 default = gradient_d8,
 options = [gradient_d8 gradient_d4],
 type = string,
-description = Method to use for calculating the slope and aspect. gradient_d8 uses 3 by 3 finite 
+description = Method to use for calculating the slope and aspect. gradient_d8 uses 3 by 3 finite
 difference window and gradient_d4 uses a two cell finite difference for x and y which mimics
 the IPW gradient function
 
+northern_hemisphere:
+default = True,
+type = bool,
+description = Boolean describing whether the model domain is in the northern hemisphere or not
+
 ################################################################################
 # Configuration for TIME section
 ################################################################################
@@ -200,15 +195,6 @@ netcdf_file:
 type =  CriticalFilename,
 description = Path to the netCDF file containing weather data
 
-# TODO pull these values from the ESPG code or it's specified in the topo section for ipw
-zone_number:
-type = int,
-description = Forces UTM zone when converting latitude and longitude to X and Y UTM coordinates
-
-# TODO pull these values from the ESPG code or it's specified in the topo section for ipw
-zone_letter:
-description = Forces UTM zone letter when converting latitude and longitude to X and Y UTM coordinates
-
 
 ################################################################################
 # air temp distribution configuration

smrf/framework/changelog.ini

@@ -6,6 +6,8 @@ date: 09-04-2019
 [changes]
 
 # Topo
+topo/basin_lon -> REMOVED
+topo/basin_lat -> REMOVED
 topo/threading -> REMOVED
 topo/dem -> REMOVED
 topo/veg_type -> REMOVED
@@ -22,8 +24,8 @@ stations -> REMOVED
 # Gridded
 gridded/file -> gridded/wrf_file
 gridded/directory -> gridded/hrrr_directory
-; gridded/zone_number -> REMOVED
-; gridded/zone_letter -> REMOVED
+gridded/zone_number -> REMOVED
+gridded/zone_letter -> REMOVED
 gridded/n_forecast_hours -> REMOVED
 
 # output

smrf/framework/model_framework.py

@@ -553,8 +553,8 @@ def distributeData_single(self):
             self._logger.info('Distributing time step %s' % t)
             # 0.1 sun angle for time step
             cosz, azimuth, rad_vec = sunang.sunang(t.astimezone(pytz.utc),
-                                            self.topo.topoConfig['basin_lat'],
-                                            self.topo.topoConfig['basin_lon'])
+                                            self.topo.basin_lat,
+                                            self.topo.basin_long)
 
             # 0.2 illumination angle
             illum_ang = None
@@ -717,8 +717,8 @@ def create_distributed_threads(self):
         t.append(Thread(target=sunang.sunang_thread,
                         name='sun_angle',
                         args=(q, self.date_time,
-                              self.topo.topoConfig['basin_lat'],
-                              self.topo.topoConfig['basin_lon'])))
+                              self.topo.basin_lat,
+                              self.topo.basin_long)))
 
         # 0.2 illumination angle
         t.append(Thread(target=radiation.shade_thread,

smrf/framework/recipes.ini

@@ -9,17 +9,7 @@ trigger:
   has_section = topo
 
 topo:
-  basin_lat = default,
-  basin_lon = default,
-  gradient_method = default
-
-[topo_filename_recipe]
-trigger_type:
-  has_value = [topo type netcdf]
-
-topo:
-  type = default,
-  filename = default
+  apply_defaults = True
 
 # TIME RECIPES
 [time_recipe]

smrf/utils/utils.py

@@ -237,7 +237,7 @@ def backup_input(data, config_obj):
         backup_config_obj.cfg['csv'][k] = fname
 
     # Copy topo files over to backup
-    ignore = ['basin_lon', 'basin_lat', 'type', 'gradient_method']
+    ignore = ['northern_hemisphere', 'gradient_method']
     for s in backup_config_obj.cfg['topo'].keys():
         src = backup_config_obj.cfg['topo'][s]
         # make not a list if lenth is 1
@@ -464,3 +464,36 @@ def grid_interpolate_deconstructed(tri, values, grid_points, method='linear'):
         return CloughTocher2DInterpolator(tri, values)(grid_points)
     elif method == 'linear':
         return LinearNDInterpolator(tri, values)(grid_points)
+
+def get_center(ds, mask_name=None):
+    '''
+    Function returns the basin center in the native coordinates of the
+    a netcdf object.
+
+    The incoming data set must contain at least and x, y and optionally
+    whatever mask name the user would like to use for calculating .
+    If no mask name is provided then the entire domain is used.
+
+    Args:
+        ds: netCDF4.Dataset object containing at least x,y, optionally
+                a mask variable name
+        mask_name: variable name in the dataset that is a mask where 1 is in
+                  the mask
+    Returns:
+        tuple: x,y of the data center in the datas native coordinates
+    '''
+    x = ds.variables['x'][:]
+    y = ds.variables['y'][:]
+
+    # Calculate the center of the basin
+    if mask_name is not None:
+        mask_id = np.argwhere(ds.variables[mask_name][:] == 1)
+
+        # Tuple is required for an upcoming deprecation in numpy
+        idx = tuple([mask_id[:, 1]])
+        idy = tuple([mask_id[:, 0]])
+
+        x = x[idx]
+        y = y[idy]
+
+    return x.mean(), y.mean()

tests/RME/config.ini

@@ -16,8 +16,6 @@
 ################################################################################
 
 [topo]
-basin_lon:                     -116.7547
-basin_lat:                     43.067
 filename:                      ./topo/topo.nc
 
 

tests/test_base_config.ini

@@ -16,8 +16,6 @@
 ################################################################################
 
 [topo]
-basin_lon:                     -116.7547
-basin_lat:                     43.067
 filename:                      ./RME/topo/topo.nc
 
 

tests/test_data_gridded.py

@@ -42,8 +42,6 @@ def test_grid_wrf(self):
             'gridded': {
                 'data_type': 'wrf',
                 'wrf_file': './RME/gridded/WRF_test.nc',
-                'zone_number': '11',
-                'zone_letter': 'N'
             },
             'system': {
                 'threading': 'False',
@@ -91,8 +89,6 @@ def test_grid_hrrr_local(self):
             'gridded': {
                 'data_type': 'hrrr_grib',
                 'hrrr_directory': './RME/gridded/hrrr_test/',
-                'zone_number': 11,
-                'zone_letter': 'N'
             },
             'time': {
                 'start_date': '2018-07-22 16:00',
@@ -152,8 +148,6 @@ def test_grid_netcdf(self):
 
         generic_grid = {'data_type': 'netcdf',
                         'netcdf_file': './RME/gridded/netcdf_test.nc',
-                        'zone_number': '11',
-                        'zone_letter': 'N',
                         'air_temp': 'air_temp',
                         'vapor_pressure': 'vapor_pressure',
                         'precip': 'precip',

tests/test_load_data.py

@@ -16,6 +16,10 @@
 from smrf.framework.model_framework import can_i_run_smrf, run_smrf
 from tests.nwrc_check import NWRCCheck
 from tests.smrf_test_case import SMRFTestCase
+import smrf
+from os.path import join, dirname, abspath
+from smrf.data import loadTopo
+import numpy as np
 
 
 @unittest.skipUnless(
@@ -146,7 +150,6 @@ def test_grid_wrf(self):
 
         adj_config = {'gridded': {'data_type': 'wrf',
                                   'wrf_file': './RME/gridded/WRF_test.nc',
-                                  'zone_number': '11',
                                   'zone_letter': 'N'},
                       'system': {'threading': 'False',
                                  'log_file': './output/log.txt'},
@@ -307,6 +310,31 @@ def test_grid_netcdf(self):
 
         self.assertIsNone(run_smrf(config))
 
+class TestLoadTopo(unittest.TestCase):
+    @classmethod
+    def setUp(self):
+        base = dirname(smrf.__file__)
+        self.test_dir = abspath(join(base, '../', 'tests'))
+        topo_config = {
+            'filename': join(self.test_dir, 'RME/topo/topo.nc'),
+            'northern_hemisphere':True,
+        }
+
+        self.topo = loadTopo.topo(
+            topo_config,
+            calcInput=False,
+            tempDir= join(self.test_dir, 'RME/output')
+        )
+
+    def test_auto_calc_lat_lon(self):
+        '''
+        Test we calculate the basin lat long correctly
+        '''
+        # Original RME
+        # basin_lon:                     -116.7547
+        # basin_lat:                     43.067
+        self.assertTrue(self.topo.basin_lat == 43.06475372378507)
+        self.assertTrue(self.topo.basin_long == -116.75395420397061)
 
 if __name__ == '__main__':
     unittest.main()