USDA-ARS-NWRC / smrf / 326

pending completion

Build # 326

Build Type

Pull #170

travis-ci-com

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/d078655acff24170422206a7dc9b585f21205ec5">d078655ac<a href="https://github.com/USDA-ARS-NWRC/smrf/commit/d078655acff24170422206a7dc9b585f21205ec5">">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/86c5a26f2cdb95bafe969cf310136c901aa64e6d">86c5a26f2</a>">86c5a26f2</a><a href="https://github.com/USDA-ARS-NWRC/smrf/commit/d078655acff24170422206a7dc9b585f21205ec5"> into 4e1ebc5df">4e1ebc5df</a>

Pull Request Pull Request #170: SMRF package

Run Details

251 of 251 new or added lines in 51 files covered. (100.0%)

3396 of 4311 relevant lines covered (78.78%)

0.79 hits per line

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

9.31

/smrf/utils/wind/model.py

from __future__ import print_function

import os
# import matplotlib.pyplot as plt
# import progressbar
from datetime import datetime

import netCDF4 as nc
import numpy as np

from smrf.utils.wind import wind_c


class wind_model():
    """

    Estimating wind speed and direction is complex terrain can be difficult due
    to the interaction of the local topography with the wind. The methods
    described here follow the work developed by Winstral and Marks (2002) and
    Winstral et al. (2009) :cite:`Winstral&Marks:2002` :cite:`Winstral&al:2009`
    which parameterizes the terrain based on the upwind direction. The
    underlying method calulates the maximum upwind slope (maxus) within a
    search distance to determine if a cell is sheltered or exposed.

    The azimuth **A** is the direction of the prevailing wind for which the
    maxus value will be calculated within a maximum search distance **dmax**.
    The maxus (**Sx**) parameter can then be estimated as the maximum value of
    the slope from the cell of interest to all of the grid cells along the
    search vector. The efficiency in selection of the maximum value can be
    increased by using the techniques from the horizon functio which calculates
    the horizon for each pixel. Therefore, less calculations can be performed.
    Negative **Sx** values indicate an exposed pixel location (shelter pixel
    was lower) and positive **Sx** values indicate a sheltered pixel (shelter
    pixel was higher).

    After all the upwind direction are calculated, the average **Sx** over a
    window is calculated. The average **Sx** accounts for larger lanscape
    obsticles that may be adjacent to the upwind direction and affect the flow.
    A window size in degrees takes the average of all **Sx**.

    Args:
        x: array of x locations
        y: array of y locations
        dem: matrix of the dem elevation values
        nthread: number of threads to use for maxus calculation

    """

    def __init__(self, x, y, dem, nthreads=1):

        self.x = x
        self.y = y
        self.dem = dem
        self.nx = len(x)
        self.ny = len(y)
        self.ngrid = self.ny * self.nx
        self.nthreads = nthreads

        self.dx = np.abs(x[1] - x[0])
        self.dy = np.abs(y[1] - y[0])

        X, Y = np.meshgrid(np.arange(0, self.nx), np.arange(0, self.ny))
        self.X = X
        self.Y = Y
        self.shape = X.shape

    def maxus(self, dmax, inc=5, inst=2, out_file='smrf_maxus.nc'):
        """
        Calculate the maxus values

        Args:
            dmax: length of outlying upwind search vector (meters)
            inc: increment between direction calculations (degrees)
            inst: Anemometer height (meters)
            out_file: NetCDF file for output results

        Returns:
            None, outputs maxus array straight to file
        """

        if (dmax % self.dx != 0):
            raise ValueError('dmax must divide evenly into the DEM')

        self.dmax = dmax
        self.inc = inc
        self.inst_hgt = inst

        # All angles that model will consider.
        swa = np.arange(0, 360, inc)
        self.directions = swa

        # initialize the output file
        self.out_file = out_file
        self.type = 'maxus'
        ex_att = {}
        ex_att['dmax'] = dmax
        # initialize output
        self.output_init(self.type, out_file, ex_att=ex_att)

        # run model over range in wind directions
        for i, angle in enumerate(swa):

            self.maxus_val = self.maxus_angle(angle, self.dmax)
            self.output(self.type, i)

    def tbreak(self, dmax, sepdist, inc=5, inst=2, out_file='smrf_tbreak.nc'):
        """
        Calculate the topobreak values

        Args:
            dmax: length of outlying upwind search vector (meters)
            sepdist: length of local max upwind slope search vector (meters)
            angle: middle upwind direction around which to run model (degrees)
            inc: increment between direction calculations (degrees)
            inst: Anemometer height (meters)
            out_file: NetCDF file for output results

        Returns:
            None, outputs maxus array straight to file

        """

        if (sepdist % self.dx != 0) | (dmax % self.dx != 0):
            raise ValueError(
                'sepdist and dmax must divide evenly into the DEM')

        self.dmax = dmax
        self.sepdist = sepdist
        self.inc = inc
        self.inst_hgt = inst

        # All angles that model will consider.
        swa = np.arange(0, 360, inc)
        self.directions = swa

        # initialize the output file
        self.out_file = out_file
        self.type = 'tbreak'
        # extra attributes
        ex_att = {}
        ex_att['dmax'] = dmax
        ex_att['sepdist'] = sepdist
        # initialize output
        self.output_init(self.type, out_file, ex_att=ex_att)

        # run model over range in wind directions
        for i, angle in enumerate(swa):

            # calculate the maxus value
            maxus_outlying = self.maxus_angle(angle, self.dmax)

            # calculate the local maxus value
            maxus_local = self.maxus_angle(angle, self.sepdist)

            self.maxus_val = maxus_local - maxus_outlying

            self.output(self.type, i)

    def maxus_angle(self, angle, dmax):
        """
        Calculate the maxus for a single direction for a search distance dmax

        Note:
            This will produce different results than the original maxus
            program. The differences are due to:

            1. Using dtype=double for the elevations
            2. Using different type of search method to find the endpoints.

            However, if the elevations are rounded to integers, the cardinal
            directions will reproduce the original results.

        Args:
            angle: middle upwind direction around which to run model (degrees)
            dmax: length of outlying upwind search vector (meters)

        Returns:
            maxus: array of maximum upwind slope values within dmax

        """

        print("Calculating maxus for direction: {}".format(angle))

        angle *= np.pi / 180

        # calculate the endpoints
        # accually use the distances to ensure that we are searching far enough
        Xi = self.X*self.dx + dmax * np.cos(angle-np.pi/2)
        Yi = self.Y*self.dy + dmax * np.sin(angle-np.pi/2)

        self.Xi = np.floor(Xi/self.dx + 0.5)
        self.Yi = np.floor(Yi/self.dy + 0.5)

        # underlying C code similar to Adams
        maxus = wind_c.call_maxus(self.x, self.y, self.dem, self.X, self.Y,
                                  self.Xi, self.Yi, self.inst_hgt,
                                  self.nthreads)

#         # my interpretation of the calculations in Python form
#         maxus = np.zeros((self.ngrid,))
#         pbar = progressbar.ProgressBar(max_value=self.ngrid)
#         j = 0
#         for index in range(5000, self.ngrid):
#             maxus[index] = self.find_maxus(index)
#             j += 1
#             pbar.update(j)
#             if j > 4999:
#                 break
#         pbar.finish()
#         maxus = np.reshape(maxus, self.shape)

        # correct for values that are their own horizon
        maxus[maxus <= -89.0] = 0

        return maxus

    def windower(self, maxus_file, window_width, wtype):
        """
        Take the maxus output and average over the window width

        Args:
            maxus_file: location of the previously calculated maxus values
            window_width: window width about the wind direction
            wtype: type of wind calculation 'maxus' or 'tbreak'

        Return:
            New file containing the windowed values
        """

        # open the previous file and get the directions
        n = nc.Dataset(maxus_file, 'r')
        directions = n.variables['direction'][:]
        self.directions = directions

        # create a new file based on the old file
        name = os.path.splitext(maxus_file)
        out_file = "%s_%iwindow.nc" % (name[0], window_width)
        self.output_init(wtype, out_file)
        self.out_file = out_file

        # determine which directions are required for each single direction
        window_width /= 2
        inc = np.mean(np.diff(directions), dtype=int)

        for i, d in enumerate(directions):

            print("Windowing direction {}".format(d))

            # determine which directions to include
            window_start = d - window_width
            window_end = d + window_width

            # to ensure that it contains the end points
            sl = np.arange(window_start, window_end+1, inc)

            # correct for edge effects
            sl[sl < 0] = sl[sl < 0] + 360
            sl[sl > 360] = sl[sl > 360] - 360

            # determine the indicies to the input file
            idx = self.ismember(directions, sl)

            # grab all the data for all directions and average
            self.maxus_val = np.mean(n.variables[wtype][idx, :], axis=0)

            # put it into the output file
            self.output(wtype, i)

        n.close()

    def ismember(self, a, b):
        bind = {}
        for i, elt in enumerate(b):
            if elt not in bind:
                bind[elt] = True
        return [bind.get(itm, False) for itm in a]

    def find_maxus(self, index):
        """
        Calculate the maxus given the start and end point

        Args:
            index: index to a point in the array

        Returns:
            maxus value for the point
        """

        start_point = np.unravel_index(index, self.shape)

        # determine the points along the endpoint line
        end_point = (self.Yi[start_point], self.Xi[start_point])
        p = self.bresenham(start_point, end_point)

        # ensure the cases where it's on the edge
        p = np.delete(p, np.where(p[:, 0] < 0)[0], axis=0)
        p = np.delete(p, np.where(p[:, 1] < 0)[0], axis=0)
        p = np.delete(p, np.where(p[:, 0] > self.ny)[0], axis=0)
        p = np.delete(p, np.where(p[:, 1] > self.nx)[0], axis=0)

        # determine the relative heights along the path
        h = self.dem[p[:, 0], p[:, 1]]  # - (self.inst_hgt + self.dem[index])

        # find the horizon for each pixel along the path
        hord = self.hord(self.x[p[:, 1]], self.y[p[:, 0]], h)

        # calculate the angle to that point
        pt = p[hord[0], :]   # point that was found for horizon
        d = np.sqrt(np.power(self.x[pt[1]] - self.x[start_point[1]], 2) +
                    np.power(self.y[pt[0]] - self.y[start_point[0]], 2))

        slope = (h[hord[0]] - (h[0] + self.inst_hgt)) / d
        maxus = np.arctan(slope) * 180 / np.pi

        return maxus

    def bresenham(self, start, end):
        """
        Python implementation of the Bresenham algorthim to find
        all the pixels that a line between start and end interscet

        Args:
            start: list of start point
            end: list of end point

        Returns:
            Array path of all points between start and end
        """
#         start = list(start)
#         end = list(end)
        path = []

        x0 = start[0]
        y0 = start[1]
        x1 = end[0]
        y1 = end[1]

        steep = abs(y1 - y0) > abs(x1 - x0)
        backward = x0 > x1

        if steep:
            x0, y0 = y0, x0
            x1, y1 = y1, x1
        if backward:
            x0, x1 = x1, x0
            y0, y1 = y1, y0

        dx = x1 - x0
        dy = abs(y1 - y0)
        error = dx / 2
        y = y0

        if y0 < y1:
            ystep = 1
        else:
            ystep = -1

        for x in range(x0, x1+1):
            if steep:
                path.append((y, x))
            else:
                path.append((x, y))

            error -= dy

            if error <= 0:
                y += ystep
                error += dx

        if backward:
            path.reverse()

        return np.array(path)

    def hord(self, x, y, z):
        '''
        Calculate the horizon pixel for all z
        This mimics the simple algorthim from Dozier 1981 but
        was adapated for use in finding the maximum upwind slope

        Works backwards from the end but looks forwards for
        the horizon

        Args:
            x: x locations for the points
            y: y locations for the points
            z: elevations for the points

        Returns:
            array of the horizon index for each point

        '''

        N = len(z)  # number of points to look at
    #     offset = 1      # offset from current point to start looking

        # preallocate the h array
        h = np.zeros(N, dtype=int)
        h[N-1] = N-1
        i = N - 2

        # work backwarks from the end for the pixels
        while i >= 0:
            h[i] = i
            j = i + 1   # looking forward
            found = False

            while not found:

                d_i = np.sqrt(np.power(x[i] - x[j], 2) +
                              np.power(y[i] - y[j], 2))
                d_h = np.sqrt(np.power(x[i] - x[h[j]], 2) +
                              np.power(y[i] - y[h[j]], 2))

                pt_i = self._slope(0, z[i]+self.inst_hgt, d_i, z[j])
                pt_h = self._slope(0, z[i]+self.inst_hgt, d_h, z[h[j]])

                if (pt_i < pt_h):
                    if (j == N-1):
                        found = True
                        h[i] = j
                    else:
                        j = h[j]
                else:
                    found = True
                    if (pt_i > pt_h):
                        h[i] = j
                    else:
                        h[i] = h[j]

            i -= 1

        return h

    def _slope(self, xi, zi, xj, zj):
        '''
        Slope between the two points
        '''

        return (zj - zi) / (xj - float(xi))

    def output_init(self, ptype, filename, ex_att=None):
        """
        Initialize a NetCDF file for outputing the maxus values or tbreak

        Args:
            ptype: type of calculation that will be saved, either 'maxus' or
                'tbreak'
            filename: filename to save the output into
            ex_att:   extra attributes to add
        """

        if ptype == 'maxus':
            var = 'maxus'
            desc = 'Maximum upwind slope'

        elif ptype == 'tbreak':
            var = 'tbreak'
            desc = 'tbreak'

        else:
            raise ValueError('''Could not determine what to output, check type
                            value (maxus or tbreak)''')

        dimensions = ('Direction', 'y', 'x')

        s = nc.Dataset(filename, 'w', 'NETCDF4')

        s.createDimension(dimensions[0], len(self.directions))
        s.createDimension(dimensions[1], self.ny)
        s.createDimension(dimensions[2], self.nx)

        # create the variables
        s.createVariable('direction', 'i', dimensions[0])
        s.createVariable('y', 'f', dimensions[1])
        s.createVariable('x', 'f', dimensions[2])
        s.createVariable(var, 'f', dimensions)

        # define some attributes
        setattr(s.variables['y'], 'units', 'meters')
        setattr(s.variables['y'], 'description', 'UTM, north south')
        setattr(s.variables['x'], 'units', 'meters')
        setattr(s.variables['x'], 'description', 'UTM, east west')
        setattr(s.variables['direction'], 'units', 'bearing')
        setattr(s.variables['direction'], 'description',
                'Wind direction from North')
        setattr(s.variables[var], 'units', 'angle')
        setattr(s.variables[var], 'description', desc)
        setattr(s, 'dateCreated', datetime.now().isoformat())

        # set attributes
        if ex_att is not None:
            for key, value in ex_att.items():
                setattr(s, key, value)

        s.variables['y'][:] = self.y
        s.variables['x'][:] = self.x

    def output(self, ptype, index):
        """
        Output the data into the out file that has previously been initialized.

        Args:
            ptype: type of calculation that will be saved, either 'maxus'
                or 'tbreak'
            index: index into the file for where to place the output
        """

        s = nc.Dataset(self.out_file, 'r+')
        s.variables['direction'][:] = self.directions
        s.variables[ptype][index, :] = self.maxus_val
        s.close()

1	from __future__ import print_function	1✔
2
3	import os	1✔
4	# import matplotlib.pyplot as plt
5	# import progressbar
6	from datetime import datetime	1✔
7
8	import netCDF4 as nc	1✔
9	import numpy as np	1✔
10
11	from smrf.utils.wind import wind_c	1✔
12
13
14	class wind_model():	1✔
15	"""
16
17	Estimating wind speed and direction is complex terrain can be difficult due
18	to the interaction of the local topography with the wind. The methods
19	described here follow the work developed by Winstral and Marks (2002) and
20	Winstral et al. (2009) :cite:`Winstral&Marks:2002` :cite:`Winstral&al:2009`
21	which parameterizes the terrain based on the upwind direction. The
22	underlying method calulates the maximum upwind slope (maxus) within a
23	search distance to determine if a cell is sheltered or exposed.
24
25	The azimuth A is the direction of the prevailing wind for which the
26	maxus value will be calculated within a maximum search distance dmax.
27	The maxus (Sx) parameter can then be estimated as the maximum value of
28	the slope from the cell of interest to all of the grid cells along the
29	search vector. The efficiency in selection of the maximum value can be
30	increased by using the techniques from the horizon functio which calculates
31	the horizon for each pixel. Therefore, less calculations can be performed.
32	Negative Sx values indicate an exposed pixel location (shelter pixel
33	was lower) and positive Sx values indicate a sheltered pixel (shelter
34	pixel was higher).
35
36	After all the upwind direction are calculated, the average Sx over a
37	window is calculated. The average Sx accounts for larger lanscape
38	obsticles that may be adjacent to the upwind direction and affect the flow.
39	A window size in degrees takes the average of all Sx.
40
41	Args:
42	x: array of x locations
43	y: array of y locations
44	dem: matrix of the dem elevation values
45	nthread: number of threads to use for maxus calculation
46
47	"""
48
49	def __init__(self, x, y, dem, nthreads=1):	1✔
50
51	self.x = x	×
52	self.y = y	×
53	self.dem = dem	×
54	self.nx = len(x)	×
55	self.ny = len(y)	×
56	self.ngrid = self.ny * self.nx	×
57	self.nthreads = nthreads	×
58
59	self.dx = np.abs(x[1] - x[0])	×
60	self.dy = np.abs(y[1] - y[0])	×
61
62	X, Y = np.meshgrid(np.arange(0, self.nx), np.arange(0, self.ny))	×
63	self.X = X	×
64	self.Y = Y	×
65	self.shape = X.shape	×
66
67	def maxus(self, dmax, inc=5, inst=2, out_file='smrf_maxus.nc'):	1✔
68	"""
69	Calculate the maxus values
70
71	Args:
72	dmax: length of outlying upwind search vector (meters)
73	inc: increment between direction calculations (degrees)
74	inst: Anemometer height (meters)
75	out_file: NetCDF file for output results
76
77	Returns:
78	None, outputs maxus array straight to file
79	"""
80
81	if (dmax % self.dx != 0):	×
82	raise ValueError('dmax must divide evenly into the DEM')	×
83
84	self.dmax = dmax	×
85	self.inc = inc	×
86	self.inst_hgt = inst	×
87
88	# All angles that model will consider.
89	swa = np.arange(0, 360, inc)	×
90	self.directions = swa	×
91
92	# initialize the output file
93	self.out_file = out_file	×
94	self.type = 'maxus'	×
95	ex_att = {}	×
96	ex_att['dmax'] = dmax	×
97	# initialize output
98	self.output_init(self.type, out_file, ex_att=ex_att)	×
99
100	# run model over range in wind directions
101	for i, angle in enumerate(swa):	×
102
103	self.maxus_val = self.maxus_angle(angle, self.dmax)	×
104	self.output(self.type, i)	×
105
106	def tbreak(self, dmax, sepdist, inc=5, inst=2, out_file='smrf_tbreak.nc'):	1✔
107	"""
108	Calculate the topobreak values
109
110	Args:
111	dmax: length of outlying upwind search vector (meters)
112	sepdist: length of local max upwind slope search vector (meters)
113	angle: middle upwind direction around which to run model (degrees)
114	inc: increment between direction calculations (degrees)
115	inst: Anemometer height (meters)
116	out_file: NetCDF file for output results
117
118	Returns:
119	None, outputs maxus array straight to file
120
121	"""
122
123	if (sepdist % self.dx != 0) \| (dmax % self.dx != 0):	×
124	raise ValueError(	×
125	'sepdist and dmax must divide evenly into the DEM')
126
127	self.dmax = dmax	×
128	self.sepdist = sepdist	×
129	self.inc = inc	×
130	self.inst_hgt = inst	×
131
132	# All angles that model will consider.
133	swa = np.arange(0, 360, inc)	×
134	self.directions = swa	×
135
136	# initialize the output file
137	self.out_file = out_file	×
138	self.type = 'tbreak'	×
139	# extra attributes
140	ex_att = {}	×
141	ex_att['dmax'] = dmax	×
142	ex_att['sepdist'] = sepdist	×
143	# initialize output
144	self.output_init(self.type, out_file, ex_att=ex_att)	×
145
146	# run model over range in wind directions
147	for i, angle in enumerate(swa):	×
148
149	# calculate the maxus value
150	maxus_outlying = self.maxus_angle(angle, self.dmax)	×
151
152	# calculate the local maxus value
153	maxus_local = self.maxus_angle(angle, self.sepdist)	×
154
155	self.maxus_val = maxus_local - maxus_outlying	×
156
157	self.output(self.type, i)	×
158
159	def maxus_angle(self, angle, dmax):	1✔
160	"""
161	Calculate the maxus for a single direction for a search distance dmax
162
163	Note:
164	This will produce different results than the original maxus
165	program. The differences are due to:
166
167	1. Using dtype=double for the elevations
168	2. Using different type of search method to find the endpoints.
169
170	However, if the elevations are rounded to integers, the cardinal
171	directions will reproduce the original results.
172
173	Args:
174	angle: middle upwind direction around which to run model (degrees)
175	dmax: length of outlying upwind search vector (meters)
176
177	Returns:
178	maxus: array of maximum upwind slope values within dmax
179
180	"""
181
182	print("Calculating maxus for direction: {}".format(angle))	×
183
184	angle *= np.pi / 180	×
185
186	# calculate the endpoints
187	# accually use the distances to ensure that we are searching far enough
188	Xi = self.Xself.dx + dmax np.cos(angle-np.pi/2)	×
189	Yi = self.Yself.dy + dmax np.sin(angle-np.pi/2)	×
190
191	self.Xi = np.floor(Xi/self.dx + 0.5)	×
192	self.Yi = np.floor(Yi/self.dy + 0.5)	×
193
194	# underlying C code similar to Adams
195	maxus = wind_c.call_maxus(self.x, self.y, self.dem, self.X, self.Y,	×
196	self.Xi, self.Yi, self.inst_hgt,
197	self.nthreads)
198
199	# # my interpretation of the calculations in Python form
200	# maxus = np.zeros((self.ngrid,))
201	# pbar = progressbar.ProgressBar(max_value=self.ngrid)
202	# j = 0
203	# for index in range(5000, self.ngrid):
204	# maxus[index] = self.find_maxus(index)
205	# j += 1
206	# pbar.update(j)
207	# if j > 4999:
208	# break
209	# pbar.finish()
210	# maxus = np.reshape(maxus, self.shape)
211
212	# correct for values that are their own horizon
213	maxus[maxus <= -89.0] = 0	×
214
215	return maxus	×
216
217	def windower(self, maxus_file, window_width, wtype):	1✔
218	"""
219	Take the maxus output and average over the window width
220
221	Args:
222	maxus_file: location of the previously calculated maxus values
223	window_width: window width about the wind direction
224	wtype: type of wind calculation 'maxus' or 'tbreak'
225
226	Return:
227	New file containing the windowed values
228	"""
229
230	# open the previous file and get the directions
231	n = nc.Dataset(maxus_file, 'r')	×
232	directions = n.variables['direction'][:]	×
233	self.directions = directions	×
234
235	# create a new file based on the old file
236	name = os.path.splitext(maxus_file)	×
237	out_file = "%s_%iwindow.nc" % (name[0], window_width)	×
238	self.output_init(wtype, out_file)	×
239	self.out_file = out_file	×
240
241	# determine which directions are required for each single direction
242	window_width /= 2	×
243	inc = np.mean(np.diff(directions), dtype=int)	×
244
245	for i, d in enumerate(directions):	×
246
247	print("Windowing direction {}".format(d))	×
248
249	# determine which directions to include
250	window_start = d - window_width	×
251	window_end = d + window_width	×
252
253	# to ensure that it contains the end points
254	sl = np.arange(window_start, window_end+1, inc)	×
255
256	# correct for edge effects
257	sl[sl < 0] = sl[sl < 0] + 360	×
258	sl[sl > 360] = sl[sl > 360] - 360	×
259
260	# determine the indicies to the input file
261	idx = self.ismember(directions, sl)	×
262
263	# grab all the data for all directions and average
264	self.maxus_val = np.mean(n.variables[wtype][idx, :], axis=0)	×
265
266	# put it into the output file
267	self.output(wtype, i)	×
268
269	n.close()	×
270
271	def ismember(self, a, b):	1✔
272	bind = {}	×
273	for i, elt in enumerate(b):	×
274	if elt not in bind:	×
275	bind[elt] = True	×
276	return [bind.get(itm, False) for itm in a]	×
277
278	def find_maxus(self, index):	1✔
279	"""
280	Calculate the maxus given the start and end point
281
282	Args:
283	index: index to a point in the array
284
285	Returns:
286	maxus value for the point
287	"""
288
289	start_point = np.unravel_index(index, self.shape)	×
290
291	# determine the points along the endpoint line
292	end_point = (self.Yi[start_point], self.Xi[start_point])	×
293	p = self.bresenham(start_point, end_point)	×
294
295	# ensure the cases where it's on the edge
296	p = np.delete(p, np.where(p[:, 0] < 0)[0], axis=0)	×
297	p = np.delete(p, np.where(p[:, 1] < 0)[0], axis=0)	×
298	p = np.delete(p, np.where(p[:, 0] > self.ny)[0], axis=0)	×
299	p = np.delete(p, np.where(p[:, 1] > self.nx)[0], axis=0)	×
300
301	# determine the relative heights along the path
302	h = self.dem[p[:, 0], p[:, 1]] # - (self.inst_hgt + self.dem[index])	×
303
304	# find the horizon for each pixel along the path
305	hord = self.hord(self.x[p[:, 1]], self.y[p[:, 0]], h)	×
306
307	# calculate the angle to that point
308	pt = p[hord[0], :] # point that was found for horizon	×
309	d = np.sqrt(np.power(self.x[pt[1]] - self.x[start_point[1]], 2) +	×
310	np.power(self.y[pt[0]] - self.y[start_point[0]], 2))
311
312	slope = (h[hord[0]] - (h[0] + self.inst_hgt)) / d	×
313	maxus = np.arctan(slope) * 180 / np.pi	×
314
315	return maxus	×
316
317	def bresenham(self, start, end):	1✔
318	"""
319	Python implementation of the Bresenham algorthim to find
320	all the pixels that a line between start and end interscet
321
322	Args:
323	start: list of start point
324	end: list of end point
325
326	Returns:
327	Array path of all points between start and end
328	"""
329	# start = list(start)
330	# end = list(end)
331	path = []	×
332
333	x0 = start[0]	×
334	y0 = start[1]	×
335	x1 = end[0]	×
336	y1 = end[1]	×
337
338	steep = abs(y1 - y0) > abs(x1 - x0)	×
339	backward = x0 > x1	×
340
341	if steep:	×
342	x0, y0 = y0, x0	×
343	x1, y1 = y1, x1	×
344	if backward:	×
345	x0, x1 = x1, x0	×
346	y0, y1 = y1, y0	×
347
348	dx = x1 - x0	×
349	dy = abs(y1 - y0)	×
350	error = dx / 2	×
351	y = y0	×
352
353	if y0 < y1:	×
354	ystep = 1	×
355	else:
356	ystep = -1	×
357
358	for x in range(x0, x1+1):	×
359	if steep:	×
360	path.append((y, x))	×
361	else:
362	path.append((x, y))	×
363
364	error -= dy	×
365
366	if error <= 0:	×
367	y += ystep	×
368	error += dx	×
369
370	if backward:	×
371	path.reverse()	×
372
373	return np.array(path)	×
374
375	def hord(self, x, y, z):	1✔
376	'''
377	Calculate the horizon pixel for all z
378	This mimics the simple algorthim from Dozier 1981 but
379	was adapated for use in finding the maximum upwind slope
380
381	Works backwards from the end but looks forwards for
382	the horizon
383
384	Args:
385	x: x locations for the points
386	y: y locations for the points
387	z: elevations for the points
388
389	Returns:
390	array of the horizon index for each point
391
392	'''
393
394	N = len(z) # number of points to look at	×
395	# offset = 1 # offset from current point to start looking
396
397	# preallocate the h array
398	h = np.zeros(N, dtype=int)	×
399	h[N-1] = N-1	×
400	i = N - 2	×
401
402	# work backwarks from the end for the pixels
403	while i >= 0:	×
404	h[i] = i	×
405	j = i + 1 # looking forward	×
406	found = False	×
407
408	while not found:	×
409
410	d_i = np.sqrt(np.power(x[i] - x[j], 2) +	×
411	np.power(y[i] - y[j], 2))
412	d_h = np.sqrt(np.power(x[i] - x[h[j]], 2) +	×
413	np.power(y[i] - y[h[j]], 2))
414
415	pt_i = self._slope(0, z[i]+self.inst_hgt, d_i, z[j])	×
416	pt_h = self._slope(0, z[i]+self.inst_hgt, d_h, z[h[j]])	×
417
418	if (pt_i < pt_h):	×
419	if (j == N-1):	×
420	found = True	×
421	h[i] = j	×
422	else:
423	j = h[j]	×
424	else:
425	found = True	×
426	if (pt_i > pt_h):	×
427	h[i] = j	×
428	else:
429	h[i] = h[j]	×
430
431	i -= 1	×
432
433	return h	×
434
435	def _slope(self, xi, zi, xj, zj):	1✔
436	'''
437	Slope between the two points
438	'''
439
440	return (zj - zi) / (xj - float(xi))	×
441
442	def output_init(self, ptype, filename, ex_att=None):	1✔
443	"""
444	Initialize a NetCDF file for outputing the maxus values or tbreak
445
446	Args:
447	ptype: type of calculation that will be saved, either 'maxus' or
448	'tbreak'
449	filename: filename to save the output into
450	ex_att: extra attributes to add
451	"""
452
453	if ptype == 'maxus':	×
454	var = 'maxus'	×
455	desc = 'Maximum upwind slope'	×
456
457	elif ptype == 'tbreak':	×
458	var = 'tbreak'	×
459	desc = 'tbreak'	×
460
461	else:
462	raise ValueError('''Could not determine what to output, check type	×
463	value (maxus or tbreak)''')
464
465	dimensions = ('Direction', 'y', 'x')	×
466
467	s = nc.Dataset(filename, 'w', 'NETCDF4')	×
468
469	s.createDimension(dimensions[0], len(self.directions))	×
470	s.createDimension(dimensions[1], self.ny)	×
471	s.createDimension(dimensions[2], self.nx)	×
472
473	# create the variables
474	s.createVariable('direction', 'i', dimensions[0])	×
475	s.createVariable('y', 'f', dimensions[1])	×
476	s.createVariable('x', 'f', dimensions[2])	×
477	s.createVariable(var, 'f', dimensions)	×
478
479	# define some attributes
480	setattr(s.variables['y'], 'units', 'meters')	×
481	setattr(s.variables['y'], 'description', 'UTM, north south')	×
482	setattr(s.variables['x'], 'units', 'meters')	×
483	setattr(s.variables['x'], 'description', 'UTM, east west')	×
484	setattr(s.variables['direction'], 'units', 'bearing')	×
485	setattr(s.variables['direction'], 'description',	×
486	'Wind direction from North')
487	setattr(s.variables[var], 'units', 'angle')	×
488	setattr(s.variables[var], 'description', desc)	×
489	setattr(s, 'dateCreated', datetime.now().isoformat())	×
490
491	# set attributes
492	if ex_att is not None:	×
493	for key, value in ex_att.items():	×
494	setattr(s, key, value)	×
495
496	s.variables['y'][:] = self.y	×
497	s.variables['x'][:] = self.x	×
498
499	def output(self, ptype, index):	1✔
500	"""
501	Output the data into the out file that has previously been initialized.
502
503	Args:
504	ptype: type of calculation that will be saved, either 'maxus'
505	or 'tbreak'
506	index: index into the file for where to place the output
507	"""
508
509	s = nc.Dataset(self.out_file, 'r+')	×
510	s.variables['direction'][:] = self.directions	×
511	s.variables[ptype][index, :] = self.maxus_val	×
512	s.close()	×

USDA-ARS-NWRC / smrf / 326

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