Auto assign basin projection attributes and basin center

smrf/data/loadGrid.py

@@ -47,39 +47,22 @@ def __init__(self, dataConfig, topo, start_date, end_date,
         self.forecast_flag = forecast_flag
         self.day_hour = day_hour
         self.n_forecast_hours = n_forecast_hours
+        self.topo = topo
 
-        # degree offset for a buffer around the model domain
+        # degree offset for a buffer around the model domain in degrees
         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']
 
-        # get the bounds of the model so that only the values inside
-        # 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]
 
-        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))
+        # get the buffer gridded data domain extents in lat long
+        self.dlat, self.dlon = self.model_domain_grid()
 
-        buff = 0.1 # buffer of bounding box in degrees
-        ur += buff
-        ll -= buff
-        self.bbox = np.append(np.flipud(ll), np.flipud(ur))
+        # This is placed long, lat on purpose as thats the HRRR class expects
+        self.bbox = np.array([self.dlon[0], self.dlat[0],
+                              self.dlon[1], self.dlat[1]])
 
         self._logger = logging.getLogger(__name__)
 
@@ -94,15 +77,24 @@ def __init__(self, dataConfig, topo, start_date, end_date,
             raise Exception('Could not resolve dataType')
 
     def model_domain_grid(self):
+        '''
+        Retrieve the bounding box for the gridded data by adding a buffer to
+        the extents of the topo domain.
 
+        Returns:
+            tuple: (dlat, dlon) Domain latitude and longitude extents
+        '''
         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'])
-        dlat[1], dlon[1] = utm.to_latlon(np.max(self.x), np.max(self.y),
-                                         int(self.dataConfig['zone_number']),
-                                         self.dataConfig['zone_letter'])
+
+        # Convert the UTM extents to lat long extents
+        ur = self.get_latlon(np.max(self.topo.x), np.max(self.topo.y))
+        ll = self.get_latlon(np.min(self.topo.x), np.min(self.topo.y))
+
+        # Put into numpy arrays for convenience later
+        dlat[0], dlon[0] = ll
+        dlat[1], dlon[1] = ur
+
         # add a buffer
         dlat[0] -= self.offset
         dlat[1] += self.offset
@@ -145,7 +137,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.topo.zone_number,
             forecast=fcast,
             forecast_flag=self.forecast_flag,
             day_hour=self.day_hour)
@@ -188,13 +180,10 @@ def load_from_hrrr(self):
         self.precip = pd.DataFrame(data['precip_int'], index=idx, columns=cols)
 
         self._logger.debug('Loading solar')
-        # solar_beam = pd.DataFrame(data['solar_beam'], index=idx, columns=cols)
-        # solar_diffuse = pd.DataFrame(data['solar_diffuse'], index=idx, columns=cols)
-        # solar = solar_beam + solar_diffuse
         solar = pd.DataFrame(data['short_wave'], index=idx, columns=cols)
         self._logger.debug('Calculating cloud factor')
         self.cloud_factor = get_hrrr_cloud(solar, self.metadata, self._logger,
-                                           self.lat, self.lon)
+                                           self.topo.basin_lat, self.topo.basin_long)
 
     def load_from_netcdf(self):
         """
@@ -221,14 +210,11 @@ def load_from_netcdf(self):
         if mlat.ndim != 2 & mlon.ndim !=2:
             [mlon, mlat] = np.meshgrid(mlon, mlat)
 
-        # get that grid cells in the model domain
-        dlat, dlon = self.model_domain_grid()
-
         # get the values that are in the modeling domain
-        ind = (mlat >= dlat[0]) & \
-            (mlat <= dlat[1]) & \
-            (mlon >= dlon[0]) & \
-            (mlon <= dlon[1])
+        ind = (mlat >= self.dlat[0]) & \
+            (mlat <= self.dlat[1]) & \
+            (mlon >= self.dlon[0]) & \
+            (mlon <= self.dlon[1])
 
         mlat = mlat[ind]
         mlon = mlon[ind]
@@ -249,7 +235,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.topo.zone_number,),
                                   axis=1)
 
         self.metadata = metadata
@@ -262,19 +248,6 @@ def load_from_netcdf(self):
             [str(tm.replace(microsecond=0)) for tm in time], tz=self.time_zone
         )
 
-        # subset the times to only those needed
-#         tzinfo = pytz.timezone(self.time_zone)
-#         time = []
-#         for t in tt:
-#             time.append(t.replace(tzinfo=tzinfo))
-#         time = np.array(time)
-
-#         time_ind = (time >= pd.to_datetime(self.start_date)) & \
-#                    (time <= pd.to_datetime(self.end_date))
-#         time = time[time_ind]
-
-#         time_idx = np.where(time_ind)[0]
-
         # GET THE DATA, ONE AT A TIME
         for v in self.variables:
 
@@ -320,24 +293,17 @@ def load_from_wrf(self):
 
         self.wrf_variables = ['GLW', 'T2', 'DWPT', 'UGRD',
                               'VGRD', 'CLDFRA', 'RAINNC']
-#         self.variables = ['thermal','air_temp','dew_point','wind_speed',
-#                             'wind_direction','cloud_factor','precip']
-
-
 
         self._logger.info('Reading data coming from WRF output: {}'.format(
             self.dataConfig['wrf_file']
             ))
         f = nc.Dataset(self.dataConfig['wrf_file'])
 
-        # DETERMINE THE MODEL DOMAIN AREA IN THE GRID
-        dlat, dlon = self.model_domain_grid()
-
         # get the values that are in the modeling domain
-        ind = (f.variables['XLAT'] >= dlat[0]) & \
-            (f.variables['XLAT'] <= dlat[1]) & \
-            (f.variables['XLONG'] >= dlon[0]) & \
-            (f.variables['XLONG'] <= dlon[1])
+        ind = (f.variables['XLAT'] >= self.dlat[0]) & \
+            (f.variables['XLAT'] <= self.dlat[1]) & \
+            (f.variables['XLONG'] >= self.dlon[0]) & \
+            (f.variables['XLONG'] <= self.dlon[1])
 
         mlat = f.variables['XLAT'][:][ind]
         mlon = f.variables['XLONG'][:][ind]
@@ -358,7 +324,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.topo.zone_number,),
                                   axis=1)
 
         self.metadata = metadata
@@ -454,6 +420,23 @@ def load_from_wrf(self):
             setattr(self, v, d.tz_convert(tz=self.time_zone))
 
 
+    def get_latlon(self, utm_x, utm_y):
+        '''
+        Convert UTM coords to Latitude and longitude
+
+        Args:
+            utm_x: UTM easting in meters in the same zone/letter as the topo
+            utm_y: UTM Northing in meters in the same zone/letter as the topo
+
+        Returns:
+            tuple: (lat,lon) latitude and longitude conversion from the UTM
+                coordinates
+        '''
+
+        lat, lon = utm.to_latlon(utm_x, utm_y, self.topo.zone_number,
+                                        northern=self.topo.northern_hemisphere)
+        return lat, lon
+
 def apply_utm(s, force_zone_number):
     """
     Calculate the utm from lat/lon for a series

smrf/data/loadTopo.py

@@ -7,10 +7,10 @@
 import numpy as np
 from netCDF4 import Dataset
 from spatialnc import ipw
+from utm import to_latlon
 
 from smrf.utils import gradient
 
-
 class topo():
     """
     Class for topo images and processing those images. Images are:
@@ -109,8 +109,62 @@ 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  = self.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 get_center(self, 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()
+
     def stoporadInput(self):
         """
         Calculate the necessary input file for stoporad
@@ -185,7 +239,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,