#198

pending completion

Build # #198

Build Type

push

coveralls-python

Committed by

web-flow

Commit Message

<a href="https://github.com/USDA-ARS-NWRC/smrf/commit/<a class=hub.com/USDA-ARS-NWRC/smrf/commit/22b64f7d7e99650ae9453a0f2722e05e26c1a476">22b64f7d7<a href="https://github.com/USDA-ARS-NWRC/smrf/commit/22b64f7d7e99650ae9453a0f2722e05e26c1a476">">Merge </a><a class="double-link" href="https://github.com/USDA-ARS-NWRC/smrf/commit/<a class="double-link" href="https://github.com/USDA-ARS-NWRC/smrf/commit/70d5a1a292419dc54fd5c5602ac712e368a44168">70d5a1a29</a>">70d5a1a29</a><a href="https://github.com/USDA-ARS-NWRC/smrf/commit/22b64f7d7e99650ae9453a0f2722e05e26c1a476"> into f8cda0533">f8cda0533</a>

Run Details

20 of 20 new or added lines in 3 files covered. (100.0%)

2743 of 3880 relevant lines covered (70.7%)

0.71 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

95.96

/smrf/data/loadGrid.py

import logging
import os

import netCDF4 as nc
import numpy as np
import pandas as pd
import pytz
import utm

from smrf.envphys import phys
from smrf.envphys.radiation import get_hrrr_cloud

try:
    from weather_forecast_retrieval import hrrr
except:
    pass



class grid():
    """
    Class for loading and storing the data, either from
    a gridded dataset in:
    - NetCDF format
    - other format

    Inputs to data() are:
    - dataConfig, from the [gridded] section
    - start_date, datetime object
    - end_date, datetime object

    """

    def __init__(self, dataConfig, topo, start_date, end_date,
                 time_zone='UTC', dataType='wrf', tempDir=None,
                 forecast_flag=False, day_hour=0, n_forecast_hours=18):

        if (tempDir is None) | (tempDir == 'WORKDIR'):
            tempDir = os.environ['WORKDIR']

        self.tempDir = tempDir
        self.dataConfig = dataConfig
        self.dataType = dataType
        self.start_date = start_date
        self.end_date = end_date
        self.time_zone = time_zone
        self.forecast_flag = forecast_flag
        self.day_hour = day_hour
        self.n_forecast_hours = n_forecast_hours

        # 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']

        # 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))

        buff = 0.1 # buffer of bounding box in degrees
        ur += buff
        ll -= buff
        self.bbox = np.append(np.flipud(ll), np.flipud(ur))

        self._logger = logging.getLogger(__name__)

        # load the data
        if dataType == 'wrf':
            self.load_from_wrf()
        elif dataType == 'netcdf':
            self.load_from_netcdf()
        elif dataType == 'hrrr_grib':
            self.load_from_hrrr()
        else:
            raise Exception('Could not resolve dataType')

    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'])
        dlat[1], dlon[1] = utm.to_latlon(np.max(self.x), np.max(self.y),
                                         int(self.dataConfig['zone_number']),
                                         self.dataConfig['zone_letter'])
        # add a buffer
        dlat[0] -= self.offset
        dlat[1] += self.offset
        dlon[0] -= self.offset
        dlon[1] += self.offset

        return dlat, dlon

    def load_from_hrrr(self):
        """
        Load the data from the High Resolution Rapid Refresh (HRRR) model
        The variables returned from the HRRR class in dataframes are

            - metadata
            - air_temp
            - relative_humidity
            - precip_int
            - cloud_factor
            - wind_u
            - wind_v

        The function will take the keys and load them into the appropriate
        objects within the `grid` class. The vapor pressure will be calculated
        from the `air_temp` and `relative_humidity`. The `wind_speed` and
        `wind_direction` will be calculated from `wind_u` and `wind_v`
        """

        self._logger.info('Reading data from from HRRR directory: {}'.format(
            self.dataConfig['hrrr_directory']
            ))

        # forecast hours for each run hour
        if not self.forecast_flag:
            fcast = [0]
        else:
            fcast = range(self.n_forecast_hours + 1)

        metadata, data = hrrr.HRRR(external_logger=self._logger).get_saved_data(
            self.start_date,
            self.end_date,
            self.bbox,
            output_dir=self.dataConfig['hrrr_directory'],
            force_zone_number=self.force_zone_number,
            forecast=fcast,
            forecast_flag=self.forecast_flag,
            day_hour=self.day_hour)

        # the data may be returned as type=object, convert to numeric
        # correct for the timezone
        for key in data.keys():
            data[key] = data[key].apply(pd.to_numeric)
            data[key] = data[key].tz_localize(tz=self.time_zone)

        self.metadata = metadata

        idx = data['air_temp'].index
        cols = data['air_temp'].columns

        self._logger.debug('Loading air_temp')
        self.air_temp = data['air_temp']

        # calculate vapor pressure
        self._logger.debug('Loading vapor_pressure')
        vp = phys.rh2vp(data['air_temp'].values, data['relative_humidity'].values)
        self.vapor_pressure = pd.DataFrame(vp, index=idx, columns=cols)

        # calculate the wind speed and wind direction
        self._logger.debug('Loading wind_speed and wind_direction')
        min_speed = 0.47

        # calculate the wind speed
        s = np.sqrt(data['wind_u']**2 + data['wind_v']**2)
        s[s < min_speed] = min_speed

        # calculate the wind direction
        d = np.degrees(np.arctan2(data['wind_v'], data['wind_u']))
        ind = d < 0
        d[ind] = d[ind] + 360
        self.wind_speed = pd.DataFrame(s, index=idx, columns=cols)
        self.wind_direction = pd.DataFrame(d, index=idx, columns=cols)

        self._logger.debug('Loading precip')
        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)

    def load_from_netcdf(self):
        """
        Load the data from a generic netcdf file

        Args:
            lat: latitude field in file, 1D array
            lon: longitude field in file, 1D array
            elev: elevation field in file, 2D array
            variable: variable name in file, 3D array
        """

        self._logger.info('Reading data coming from netcdf: {}'.format(
                            self.dataConfig['netcdf_file'])
                          )

        f = nc.Dataset(self.dataConfig['netcdf_file'], 'r')

        # GET THE LAT, LON, ELEV FROM THE FILE
        mlat = f.variables['lat'][:]
        mlon = f.variables['lon'][:]
        mhgt = f.variables['elev'][:]

        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])

        mlat = mlat[ind]
        mlon = mlon[ind]
        mhgt = mhgt[ind]

        # GET THE METADATA
        # create some fake station names based on the index
        a = np.argwhere(ind)
        primary_id = ['grid_y%i_x%i' % (i[0], i[1]) for i in a]
        self._logger.debug('{} grid cells within model domain'.format(len(a)))

        # create a metadata dataframe to store all the grid info
        metadata = pd.DataFrame(index=primary_id,
                                columns=('X', 'Y', 'latitude',
                                         'longitude', 'elevation'))

        metadata['latitude'] = mlat.flatten()
        metadata['longitude'] = mlon.flatten()
        metadata['elevation'] = mhgt.flatten()
        metadata = metadata.apply(apply_utm,
                                  args=(self.force_zone_number,),
                                  axis=1)

        self.metadata = metadata

        # GET THE TIMES
        t = f.variables['time']
        time = nc.num2date(t[:].astype(int), t.getncattr('units'), t.getncattr('calendar'))
        # Drop milliseconds and prepare to use as pandas DataFrame index
        time = pd.DatetimeIndex(
            [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:

            if v in self.dataConfig:
                v_file = self.dataConfig[v]
                self._logger.debug('Loading {} from {}'.format(v, v_file))

                df = pd.DataFrame(index=time, columns=primary_id)
                for i in a:
                    g = 'grid_y%i_x%i' % (i[0], i[1])
                    df[g] = f.variables[v_file][:, i[0], i[1]]

                # deal with any fillValues
                try:
                    fv = f.variables[v_file].getncattr('_FillValue')
                    df.replace(fv, np.nan, inplace=True)
                except:
                    pass

                # Set variable and subset by start and end time
                setattr(self, v, df[self.start_date:self.end_date])

    def load_from_wrf(self):
        """
        Load the data from a netcdf file. This was setup to work with a WRF
        output file, i.e. wrf_out so it's going to look for the following
        variables:
        - Times
        - XLAT
        - XLONG
        - HGT
        - T2
        - DWPT
        - GLW
        - RAINNC
        - CLDFRA
        - UGRD
        - VGRD

        Each cell will be identified by grid_IX_IY

        """

        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])

        mlat = f.variables['XLAT'][:][ind]
        mlon = f.variables['XLONG'][:][ind]
        mhgt = f.variables['HGT'][:][ind]

        # GET THE METADATA
        # create some fake station names based on the index
        a = np.argwhere(ind)
        primary_id = ['grid_y%i_x%i' % (i[0], i[1]) for i in a]
        self._logger.debug('{} grid cells within model domain'.format(len(a)))

        # create a metadata dataframe to store all the grid info
        metadata = pd.DataFrame(index=primary_id,
                                columns=('X', 'Y', 'latitude',
                                         'longitude', 'elevation'))

        metadata['latitude'] = mlat.flatten()
        metadata['longitude'] = mlon.flatten()
        metadata['elevation'] = mhgt.flatten()
        metadata = metadata.apply(apply_utm,
                                  args=(self.force_zone_number,),
                                  axis=1)

        self.metadata = metadata

        # GET THE TIMES
        t = f.variables['Times']
        t.set_auto_maskandscale(False)
        try:
            times = [('').join(v) for v in t]
        except TypeError:
            times = []
            for v in t:
                times.append(''.join([s.decode('utf-8') for s in v]))
        except Exception:
            raise Exception('Could not convert WRF times to readable format')

        times = [v.replace('_', ' ') for v in times]  # remove the underscore
        in_utc = str(self.time_zone).lower() == str(pytz.UTC).lower()
        time = pd.to_datetime(times, utc=in_utc)

        # subset the times to only those needed
        time_ind = (time >= self.start_date) & (time <= self.end_date)
        time = time[time_ind]

        # GET THE DATA, ONE AT A TIME
        self._logger.debug('Loading air_temp')
        self.air_temp = pd.DataFrame(index=time, columns=primary_id)
        for i in a:
            g = 'grid_y%i_x%i' % (i[0], i[1])
            v = f.variables['T2'][time_ind, i[0], i[1]] - 273.15
            self.air_temp[g] = v

        self._logger.debug('Loading dew_point and calculating vapor_pressure')
        self.dew_point = pd.DataFrame(index=time, columns=primary_id)
        for i in a:
            g = 'grid_y%i_x%i' % (i[0], i[1])
            v = f.variables['DWPT'][time_ind, i[0], i[1]]
            self.dew_point[g] = v

        self._logger.debug('Calculating vapor_pressure')
        satvp = phys.satvp(self.dew_point.values)
        self.vapor_pressure = pd.DataFrame(satvp, index=time, columns=primary_id)

        self._logger.debug('Loading thermal')
        self.thermal = pd.DataFrame(index=time, columns=primary_id)
        for i in a:
            g = 'grid_y%i_x%i' % (i[0], i[1])
            v = f.variables['GLW'][time_ind, i[0], i[1]]
            self.thermal[g] = v

        self._logger.debug('Loading cloud_factor')
        self.cloud_factor = pd.DataFrame(index=time, columns=primary_id)
        cf = 1 - np.mean(f.variables['CLDFRA'][time_ind, :], axis=1)
        for i in a:
            g = 'grid_y%i_x%i' % (i[0], i[1])
            v = cf[:, i[0], i[1]]
            self.cloud_factor[g] = v

        self._logger.debug('Loading wind_speed and wind_direction')
        self.wind_speed = pd.DataFrame(index=time, columns=primary_id)
        self.wind_direction = pd.DataFrame(index=time, columns=primary_id)
        min_speed = 0.47

        u10 = f.variables['UGRD'][time_ind, :]
        v10 = f.variables['VGRD'][time_ind, :]

        # calculate the wind speed
        s = np.sqrt(u10**2 + v10**2)
        s[s < min_speed] = min_speed

        # calculate the wind direction
        d = np.degrees(np.arctan2(v10, u10))

        ind = d < 0
        d[ind] = d[ind] + 360

        for i in a:
            g = 'grid_y%i_x%i' % (i[0], i[1])

            self.wind_speed[g] = s[:, i[0], i[1]]
            self.wind_direction[g] = d[:, i[0], i[1]]

        self._logger.debug('Loading precip')
        self.precip = pd.DataFrame(index=time, columns=primary_id)
        precip = np.diff(f.variables['RAINNC'][time_ind, :], axis=0)
        for i in a:
            g = 'grid_y%i_x%i' % (i[0], i[1])
            self.precip[g] = np.concatenate(([0], precip[:, i[0], i[1]]))

        # correct for the timezone and get only the desired dates
        for v in self.variables:
            d = getattr(self, v)
            setattr(self, v, d.tz_convert(tz=self.time_zone))


def apply_utm(s, force_zone_number):
    """
    Calculate the utm from lat/lon for a series

    Args:
        s: pandas series with fields latitude and longitude
        force_zone_number: default None, zone number to force to

    Returns:
        s: pandas series with fields 'X' and 'Y' filled
    """
    p = utm.from_latlon(s.latitude, s.longitude,
                        force_zone_number=force_zone_number)
    s['X'] = p[0]
    s['Y'] = p[1]
    return s

1	import logging	1✔
2	import os	1✔
3
4	import netCDF4 as nc	1✔
5	import numpy as np	1✔
6	import pandas as pd	1✔
7	import pytz	1✔
8	import utm	1✔
9
10	from smrf.envphys import phys	1✔
11	from smrf.envphys.radiation import get_hrrr_cloud	1✔
12
13	try:	1✔
14	from weather_forecast_retrieval import hrrr	1✔
15	except:	×
16	pass	×
17
18
19
20	class grid():	1✔
21	"""
22	Class for loading and storing the data, either from
23	a gridded dataset in:
24	- NetCDF format
25	- other format
26
27	Inputs to data() are:
28	- dataConfig, from the [gridded] section
29	- start_date, datetime object
30	- end_date, datetime object
31
32	"""
33
34	def __init__(self, dataConfig, topo, start_date, end_date,	1✔
35	time_zone='UTC', dataType='wrf', tempDir=None,
36	forecast_flag=False, day_hour=0, n_forecast_hours=18):
37
38	if (tempDir is None) \| (tempDir == 'WORKDIR'):	1✔
39	tempDir = os.environ['WORKDIR']	×
40
41	self.tempDir = tempDir	1✔
42	self.dataConfig = dataConfig	1✔
43	self.dataType = dataType	1✔
44	self.start_date = start_date	1✔
45	self.end_date = end_date	1✔
46	self.time_zone = time_zone	1✔
47	self.forecast_flag = forecast_flag	1✔
48	self.day_hour = day_hour	1✔
49	self.n_forecast_hours = n_forecast_hours	1✔
50
51	# degree offset for a buffer around the model domain
52	self.offset = 0.1	1✔
53
54	self.force_zone_number = None	1✔
55	if 'zone_number' in dataConfig:	1✔
56	self.force_zone_number = dataConfig['zone_number']	1✔
57
58	# The data that will be output
59	self.variables = ['air_temp', 'vapor_pressure', 'precip', 'wind_speed',	1✔
60	'wind_direction', 'cloud_factor', 'thermal']
61
62	# get the bounds of the model so that only the values inside
63	# the model domain are used
64	self.x = topo.x	1✔
65	self.y = topo.y	1✔
66	self.lat = topo.topoConfig['basin_lat']	1✔
67	self.lon = topo.topoConfig['basin_lon']	1✔
68
69	# get the zone number and the bounding box
70	u = utm.from_latlon(topo.topoConfig['basin_lat'],	1✔
71	topo.topoConfig['basin_lon'],
72	self.force_zone_number)
73	self.zone_number = u[2]	1✔
74	self.zone_letter = u[3]	1✔
75
76	ur = np.array(utm.to_latlon(np.max(self.x), np.max(self.y), self.zone_number, self.zone_letter))	1✔
77	ll = np.array(utm.to_latlon(np.min(self.x), np.min(self.y), self.zone_number, self.zone_letter))	1✔
78
79	buff = 0.1 # buffer of bounding box in degrees	1✔
80	ur += buff	1✔
81	ll -= buff	1✔
82	self.bbox = np.append(np.flipud(ll), np.flipud(ur))	1✔
83
84	self._logger = logging.getLogger(__name__)	1✔
85
86	# load the data
87	if dataType == 'wrf':	1✔
88	self.load_from_wrf()	1✔
89	elif dataType == 'netcdf':	1✔
90	self.load_from_netcdf()	1✔
91	elif dataType == 'hrrr_grib':	1✔
92	self.load_from_hrrr()	1✔
93	else:
94	raise Exception('Could not resolve dataType')	×
95
96	def model_domain_grid(self):	1✔
97
98	dlat = np.zeros((2,))	1✔
99	dlon = np.zeros_like(dlat)	1✔
100	dlat[0], dlon[0] = utm.to_latlon(np.min(self.x), np.min(self.y),	1✔
101	int(self.dataConfig['zone_number']),
102	self.dataConfig['zone_letter'])
103	dlat[1], dlon[1] = utm.to_latlon(np.max(self.x), np.max(self.y),	1✔
104	int(self.dataConfig['zone_number']),
105	self.dataConfig['zone_letter'])
106	# add a buffer
107	dlat[0] -= self.offset	1✔
108	dlat[1] += self.offset	1✔
109	dlon[0] -= self.offset	1✔
110	dlon[1] += self.offset	1✔
111
112	return dlat, dlon	1✔
113
114	def load_from_hrrr(self):	1✔
115	"""
116	Load the data from the High Resolution Rapid Refresh (HRRR) model
117	The variables returned from the HRRR class in dataframes are
118
119	- metadata
120	- air_temp
121	- relative_humidity
122	- precip_int
123	- cloud_factor
124	- wind_u
125	- wind_v
126
127	The function will take the keys and load them into the appropriate
128	objects within the `grid` class. The vapor pressure will be calculated
129	from the `air_temp` and `relative_humidity`. The `wind_speed` and
130	`wind_direction` will be calculated from `wind_u` and `wind_v`
131	"""
132
133	self._logger.info('Reading data from from HRRR directory: {}'.format(	1✔
134	self.dataConfig['hrrr_directory']
135	))
136
137	# forecast hours for each run hour
138	if not self.forecast_flag:	1✔
139	fcast = [0]	1✔
140	else:
141	fcast = range(self.n_forecast_hours + 1)	×
142
143	metadata, data = hrrr.HRRR(external_logger=self._logger).get_saved_data(	1✔
144	self.start_date,
145	self.end_date,
146	self.bbox,
147	output_dir=self.dataConfig['hrrr_directory'],
148	force_zone_number=self.force_zone_number,
149	forecast=fcast,
150	forecast_flag=self.forecast_flag,
151	day_hour=self.day_hour)
152
153	# the data may be returned as type=object, convert to numeric
154	# correct for the timezone
155	for key in data.keys():	1✔
156	data[key] = data[key].apply(pd.to_numeric)	1✔
157	data[key] = data[key].tz_localize(tz=self.time_zone)	1✔
158
159	self.metadata = metadata	1✔
160
161	idx = data['air_temp'].index	1✔
162	cols = data['air_temp'].columns	1✔
163
164	self._logger.debug('Loading air_temp')	1✔
165	self.air_temp = data['air_temp']	1✔
166
167	# calculate vapor pressure
168	self._logger.debug('Loading vapor_pressure')	1✔
169	vp = phys.rh2vp(data['air_temp'].values, data['relative_humidity'].values)	1✔
170	self.vapor_pressure = pd.DataFrame(vp, index=idx, columns=cols)	1✔
171
172	# calculate the wind speed and wind direction
173	self._logger.debug('Loading wind_speed and wind_direction')	1✔
174	min_speed = 0.47	1✔
175
176	# calculate the wind speed
177	s = np.sqrt(data['wind_u']2 + data['wind_v']2)	1✔
178	s[s < min_speed] = min_speed	1✔
179
180	# calculate the wind direction
181	d = np.degrees(np.arctan2(data['wind_v'], data['wind_u']))	1✔
182	ind = d < 0	1✔
183	d[ind] = d[ind] + 360	1✔
184	self.wind_speed = pd.DataFrame(s, index=idx, columns=cols)	1✔
185	self.wind_direction = pd.DataFrame(d, index=idx, columns=cols)	1✔
186
187	self._logger.debug('Loading precip')	1✔
188	self.precip = pd.DataFrame(data['precip_int'], index=idx, columns=cols)	1✔
189
190	self._logger.debug('Loading solar')	1✔
191	# solar_beam = pd.DataFrame(data['solar_beam'], index=idx, columns=cols)
192	# solar_diffuse = pd.DataFrame(data['solar_diffuse'], index=idx, columns=cols)
193	# solar = solar_beam + solar_diffuse
194	solar = pd.DataFrame(data['short_wave'], index=idx, columns=cols)	1✔
195	self._logger.debug('Calculating cloud factor')	1✔
196	self.cloud_factor = get_hrrr_cloud(solar, self.metadata, self._logger,	1✔
197	self.lat, self.lon)
198
199	def load_from_netcdf(self):	1✔
200	"""
201	Load the data from a generic netcdf file
202
203	Args:
204	lat: latitude field in file, 1D array
205	lon: longitude field in file, 1D array
206	elev: elevation field in file, 2D array
207	variable: variable name in file, 3D array
208	"""
209
210	self._logger.info('Reading data coming from netcdf: {}'.format(	1✔
211	self.dataConfig['netcdf_file'])
212	)
213
214	f = nc.Dataset(self.dataConfig['netcdf_file'], 'r')	1✔
215
216	# GET THE LAT, LON, ELEV FROM THE FILE
217	mlat = f.variables['lat'][:]	1✔
218	mlon = f.variables['lon'][:]	1✔
219	mhgt = f.variables['elev'][:]	1✔
220
221	if mlat.ndim != 2 & mlon.ndim !=2:	1✔
222	[mlon, mlat] = np.meshgrid(mlon, mlat)	×
223
224	# get that grid cells in the model domain
225	dlat, dlon = self.model_domain_grid()	1✔
226
227	# get the values that are in the modeling domain
228	ind = (mlat >= dlat[0]) & \	1✔
229	(mlat <= dlat[1]) & \
230	(mlon >= dlon[0]) & \
231	(mlon <= dlon[1])
232
233	mlat = mlat[ind]	1✔
234	mlon = mlon[ind]	1✔
235	mhgt = mhgt[ind]	1✔
236
237	# GET THE METADATA
238	# create some fake station names based on the index
239	a = np.argwhere(ind)	1✔
240	primary_id = ['grid_y%i_x%i' % (i[0], i[1]) for i in a]	1✔
241	self._logger.debug('{} grid cells within model domain'.format(len(a)))	1✔
242
243	# create a metadata dataframe to store all the grid info
244	metadata = pd.DataFrame(index=primary_id,	1✔
245	columns=('X', 'Y', 'latitude',
246	'longitude', 'elevation'))
247
248	metadata['latitude'] = mlat.flatten()	1✔
249	metadata['longitude'] = mlon.flatten()	1✔
250	metadata['elevation'] = mhgt.flatten()	1✔
251	metadata = metadata.apply(apply_utm,	1✔
252	args=(self.force_zone_number,),
253	axis=1)
254
255	self.metadata = metadata	1✔
256
257	# GET THE TIMES
258	t = f.variables['time']	1✔
259	time = nc.num2date(t[:].astype(int), t.getncattr('units'), t.getncattr('calendar'))	1✔
260	# Drop milliseconds and prepare to use as pandas DataFrame index
261	time = pd.DatetimeIndex(	1✔
262	[str(tm.replace(microsecond=0)) for tm in time], tz=self.time_zone
263	)
264
265	# subset the times to only those needed
266	# tzinfo = pytz.timezone(self.time_zone)
267	# time = []
268	# for t in tt:
269	# time.append(t.replace(tzinfo=tzinfo))
270	# time = np.array(time)
271
272	# time_ind = (time >= pd.to_datetime(self.start_date)) & \
273	# (time <= pd.to_datetime(self.end_date))
274	# time = time[time_ind]
275
276	# time_idx = np.where(time_ind)[0]
277
278	# GET THE DATA, ONE AT A TIME
279	for v in self.variables:	1✔
280
281	if v in self.dataConfig:	1✔
282	v_file = self.dataConfig[v]	1✔
283	self._logger.debug('Loading {} from {}'.format(v, v_file))	1✔
284
285	df = pd.DataFrame(index=time, columns=primary_id)	1✔
286	for i in a:	1✔
287	g = 'grid_y%i_x%i' % (i[0], i[1])	1✔
288	df[g] = f.variables[v_file][:, i[0], i[1]]	1✔
289
290	# deal with any fillValues
291	try:	1✔
292	fv = f.variables[v_file].getncattr('_FillValue')	1✔
293	df.replace(fv, np.nan, inplace=True)	×
294	except:	1✔
295	pass	1✔
296
297	# Set variable and subset by start and end time
298	setattr(self, v, df[self.start_date:self.end_date])	1✔
299
300	def load_from_wrf(self):	1✔
301	"""
302	Load the data from a netcdf file. This was setup to work with a WRF
303	output file, i.e. wrf_out so it's going to look for the following
304	variables:
305	- Times
306	- XLAT
307	- XLONG
308	- HGT
309	- T2
310	- DWPT
311	- GLW
312	- RAINNC
313	- CLDFRA
314	- UGRD
315	- VGRD
316
317	Each cell will be identified by grid_IX_IY
318
319	"""
320
321	self.wrf_variables = ['GLW', 'T2', 'DWPT', 'UGRD',	1✔
322	'VGRD', 'CLDFRA', 'RAINNC']
323	# self.variables = ['thermal','air_temp','dew_point','wind_speed',
324	# 'wind_direction','cloud_factor','precip']
325
326
327
328	self._logger.info('Reading data coming from WRF output: {}'.format(	1✔
329	self.dataConfig['wrf_file']
330	))
331	f = nc.Dataset(self.dataConfig['wrf_file'])	1✔
332
333	# DETERMINE THE MODEL DOMAIN AREA IN THE GRID
334	dlat, dlon = self.model_domain_grid()	1✔
335
336	# get the values that are in the modeling domain
337	ind = (f.variables['XLAT'] >= dlat[0]) & \	1✔
338	(f.variables['XLAT'] <= dlat[1]) & \
339	(f.variables['XLONG'] >= dlon[0]) & \
340	(f.variables['XLONG'] <= dlon[1])
341
342	mlat = f.variables['XLAT'][:][ind]	1✔
343	mlon = f.variables['XLONG'][:][ind]	1✔
344	mhgt = f.variables['HGT'][:][ind]	1✔
345
346	# GET THE METADATA
347	# create some fake station names based on the index
348	a = np.argwhere(ind)	1✔
349	primary_id = ['grid_y%i_x%i' % (i[0], i[1]) for i in a]	1✔
350	self._logger.debug('{} grid cells within model domain'.format(len(a)))	1✔
351
352	# create a metadata dataframe to store all the grid info
353	metadata = pd.DataFrame(index=primary_id,	1✔
354	columns=('X', 'Y', 'latitude',
355	'longitude', 'elevation'))
356
357	metadata['latitude'] = mlat.flatten()	1✔
358	metadata['longitude'] = mlon.flatten()	1✔
359	metadata['elevation'] = mhgt.flatten()	1✔
360	metadata = metadata.apply(apply_utm,	1✔
361	args=(self.force_zone_number,),
362	axis=1)
363
364	self.metadata = metadata	1✔
365
366	# GET THE TIMES
367	t = f.variables['Times']	1✔
368	t.set_auto_maskandscale(False)	1✔
369	try:	1✔
370	times = [('').join(v) for v in t]	1✔
371	except TypeError:	1✔
372	times = []	1✔
373	for v in t:	1✔
374	times.append(''.join([s.decode('utf-8') for s in v]))	1✔
375	except Exception:	×
376	raise Exception('Could not convert WRF times to readable format')	×
377
378	times = [v.replace('_', ' ') for v in times] # remove the underscore	1✔
379	in_utc = str(self.time_zone).lower() == str(pytz.UTC).lower()	1✔
380	time = pd.to_datetime(times, utc=in_utc)	1✔
381
382	# subset the times to only those needed
383	time_ind = (time >= self.start_date) & (time <= self.end_date)	1✔
384	time = time[time_ind]	1✔
385
386	# GET THE DATA, ONE AT A TIME
387	self._logger.debug('Loading air_temp')	1✔
388	self.air_temp = pd.DataFrame(index=time, columns=primary_id)	1✔
389	for i in a:	1✔
390	g = 'grid_y%i_x%i' % (i[0], i[1])	1✔
391	v = f.variables['T2'][time_ind, i[0], i[1]] - 273.15	1✔
392	self.air_temp[g] = v	1✔
393
394	self._logger.debug('Loading dew_point and calculating vapor_pressure')	1✔
395	self.dew_point = pd.DataFrame(index=time, columns=primary_id)	1✔
396	for i in a:	1✔
397	g = 'grid_y%i_x%i' % (i[0], i[1])	1✔
398	v = f.variables['DWPT'][time_ind, i[0], i[1]]	1✔
399	self.dew_point[g] = v	1✔
400
401	self._logger.debug('Calculating vapor_pressure')	1✔
402	satvp = phys.satvp(self.dew_point.values)	1✔
403	self.vapor_pressure = pd.DataFrame(satvp, index=time, columns=primary_id)	1✔
404
405	self._logger.debug('Loading thermal')	1✔
406	self.thermal = pd.DataFrame(index=time, columns=primary_id)	1✔
407	for i in a:	1✔
408	g = 'grid_y%i_x%i' % (i[0], i[1])	1✔
409	v = f.variables['GLW'][time_ind, i[0], i[1]]	1✔
410	self.thermal[g] = v	1✔
411
412	self._logger.debug('Loading cloud_factor')	1✔
413	self.cloud_factor = pd.DataFrame(index=time, columns=primary_id)	1✔
414	cf = 1 - np.mean(f.variables['CLDFRA'][time_ind, :], axis=1)	1✔
415	for i in a:	1✔
416	g = 'grid_y%i_x%i' % (i[0], i[1])	1✔
417	v = cf[:, i[0], i[1]]	1✔
418	self.cloud_factor[g] = v	1✔
419
420	self._logger.debug('Loading wind_speed and wind_direction')	1✔
421	self.wind_speed = pd.DataFrame(index=time, columns=primary_id)	1✔
422	self.wind_direction = pd.DataFrame(index=time, columns=primary_id)	1✔
423	min_speed = 0.47	1✔
424
425	u10 = f.variables['UGRD'][time_ind, :]	1✔
426	v10 = f.variables['VGRD'][time_ind, :]	1✔
427
428	# calculate the wind speed
429	s = np.sqrt(u102 + v102)	1✔
430	s[s < min_speed] = min_speed	1✔
431
432	# calculate the wind direction
433	d = np.degrees(np.arctan2(v10, u10))	1✔
434
435	ind = d < 0	1✔
436	d[ind] = d[ind] + 360	1✔
437
438	for i in a:	1✔
439	g = 'grid_y%i_x%i' % (i[0], i[1])	1✔
440
441	self.wind_speed[g] = s[:, i[0], i[1]]	1✔
442	self.wind_direction[g] = d[:, i[0], i[1]]	1✔
443
444	self._logger.debug('Loading precip')	1✔
445	self.precip = pd.DataFrame(index=time, columns=primary_id)	1✔
446	precip = np.diff(f.variables['RAINNC'][time_ind, :], axis=0)	1✔
447	for i in a:	1✔
448	g = 'grid_y%i_x%i' % (i[0], i[1])	1✔
449	self.precip[g] = np.concatenate(([0], precip[:, i[0], i[1]]))	1✔
450
451	# correct for the timezone and get only the desired dates
452	for v in self.variables:	1✔
453	d = getattr(self, v)	1✔
454	setattr(self, v, d.tz_convert(tz=self.time_zone))	1✔
455
456
457	def apply_utm(s, force_zone_number):	1✔
458	"""
459	Calculate the utm from lat/lon for a series
460
461	Args:
462	s: pandas series with fields latitude and longitude
463	force_zone_number: default None, zone number to force to
464
465	Returns:
466	s: pandas series with fields 'X' and 'Y' filled
467	"""
468	p = utm.from_latlon(s.latitude, s.longitude,	1✔
469	force_zone_number=force_zone_number)
470	s['X'] = p[0]	1✔
471	s['Y'] = p[1]	1✔
472	return s	1✔

USDA-ARS-NWRC / smrf / #198

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous