Skip to content
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
46 changes: 32 additions & 14 deletions rocketpy/Environment.py
Original file line number Diff line number Diff line change
Expand Up @@ -25,12 +25,27 @@

try:
import netCDF4
from netCDF4 import Dataset
except ImportError:
has_netCDF4 = False
warnings.warn(
"Unable to load netCDF4. NetCDF files and OPeNDAP will not be imported.",
ImportWarning,
)
else:
has_netCDF4 = True


def requires_netCDF4(func):
def wrapped_func(*args, **kwargs):
if has_netCDF4:
func(*args, **kwargs)
else:
raise ImportError(
"This feature requires netCDF4 to be installed. Install it with `pip install netCDF4`"
)

return wrapped_func


from .Function import Function

Expand Down Expand Up @@ -348,7 +363,7 @@ def __init__(
self.date = None

# Initialize constants
self.earthRadius = 6.3781 * (10 ** 6)
self.earthRadius = 6.3781 * (10**6)
self.airGasConstant = 287.05287 # in J/K/Kg

# Initialize atmosphere
Expand Down Expand Up @@ -487,6 +502,7 @@ def setElevation(self, elevation="Open-Elevation"):
" Open-Elevation API. See Environment.setLocation."
)

@requires_netCDF4
def setTopographicProfile(self, type, file, dictionary="netCDF4", crs=None):
"""[UNDER CONSTRUCTION] Defines the Topographic profile, importing data
from previous downloaded files. Mainly data from the Shuttle Radar
Expand Down Expand Up @@ -515,7 +531,7 @@ def setTopographicProfile(self, type, file, dictionary="netCDF4", crs=None):

if type == "NASADEM_HGT":
if dictionary == "netCDF4":
rootgrp = Dataset(file, "r", format="NETCDF4")
rootgrp = netCDF4.Dataset(file, "r", format="NETCDF4")
self.elevLonArray = rootgrp.variables["lon"][:].tolist()
self.elevLatArray = rootgrp.variables["lat"][:].tolist()
self.elevArray = rootgrp.variables["NASADEM_HGT"][:].tolist()
Expand Down Expand Up @@ -1686,6 +1702,7 @@ def processNOAARUCSounding(self, file):
# Save maximum expected height
self.maxExpectedHeight = pressure_array[-1, 0]

@requires_netCDF4
def processForecastReanalysis(self, file, dictionary):
"""Import and process atmospheric data from weather forecasts
and reanalysis given as netCDF or OPeNDAP files.
Expand Down Expand Up @@ -1944,7 +1961,7 @@ def processForecastReanalysis(self, file, dictionary):
windV = ((y2 - y) / (y2 - y1)) * f_x_y1 + ((y - y1) / (y2 - y1)) * f_x_y2

# Determine wind speed, heading and direction
windSpeed = np.sqrt(windU ** 2 + windV ** 2)
windSpeed = np.sqrt(windU**2 + windV**2)
windHeading = np.arctan2(windU, windV) * (180 / np.pi) % 360
windDirection = (windHeading - 180) % 360

Expand Down Expand Up @@ -2079,6 +2096,7 @@ def processForecastReanalysis(self, file, dictionary):

return None

@requires_netCDF4
def processEnsemble(self, file, dictionary):
"""Import and process atmospheric data from weather ensembles
given as netCDF or OPeNDAP files.
Expand Down Expand Up @@ -2361,7 +2379,7 @@ def processEnsemble(self, file, dictionary):
windV = ((y2 - y) / (y2 - y1)) * f_x_y1 + ((y - y1) / (y2 - y1)) * f_x_y2

# Determine wind speed, heading and direction
windSpeed = np.sqrt(windU ** 2 + windV ** 2)
windSpeed = np.sqrt(windU**2 + windV**2)
windHeading = np.arctan2(windU, windV) * (180 / np.pi) % 360
windDirection = (windHeading - 180) % 360

Expand Down Expand Up @@ -3249,7 +3267,7 @@ def geodesicToUtm(self, lat, lon, datum):

# Evaluate reference parameters
K0 = 1 - 1 / 2500
e2 = 2 * flattening - flattening ** 2
e2 = 2 * flattening - flattening**2
e2lin = e2 / (1 - e2)

# Evaluate auxiliary parameters
Expand All @@ -3272,9 +3290,9 @@ def geodesicToUtm(self, lat, lon, datum):

# Evaluate new auxiliary parameters
J = (1 - t + c) * ag * ag * ag / 6
K = (5 - 18 * t + t * t + 72 * c - 58 * e2lin) * (ag ** 5) / 120
K = (5 - 18 * t + t * t + 72 * c - 58 * e2lin) * (ag**5) / 120
L = (5 - t + 9 * c + 4 * c * c) * ag * ag * ag * ag / 24
M = (61 - 58 * t + t * t + 600 * c - 330 * e2lin) * (ag ** 6) / 720
M = (61 - 58 * t + t * t + 600 * c - 330 * e2lin) * (ag**6) / 720

# Evaluate the final coordinates
x = 500000 + K0 * n * (ag + J + K)
Expand Down Expand Up @@ -3347,7 +3365,7 @@ def utmToGeodesic(self, x, y, utmZone, hemis, datum):

# Calculate reference values
K0 = 1 - 1 / 2500
e2 = 2 * flattening - flattening ** 2
e2 = 2 * flattening - flattening**2
e2lin = e2 / (1 - e2)
e1 = (1 - (1 - e2) ** 0.5) / (1 + (1 - e2) ** 0.5)

Expand All @@ -3371,20 +3389,20 @@ def utmToGeodesic(self, x, y, utmZone, hemis, datum):
t1 = np.tan(lat1) ** 2
n1 = semiMajorAxis / ((1 - e2 * (np.sin(lat1) ** 2)) ** 0.5)
quoc = (1 - e2 * np.sin(lat1) * np.sin(lat1)) ** 3
r1 = semiMajorAxis * (1 - e2) / (quoc ** 0.5)
r1 = semiMajorAxis * (1 - e2) / (quoc**0.5)
d = (x - 500000) / (n1 * K0)

# Calculate other auxiliary values
I = (5 + 3 * t1 + 10 * c1 - 4 * c1 * c1 - 9 * e2lin) * d * d * d * d / 24
J = (
(61 + 90 * t1 + 298 * c1 + 45 * t1 * t1 - 252 * e2lin - 3 * c1 * c1)
* (d ** 6)
* (d**6)
/ 720
)
K = d - (1 + 2 * t1 + c1) * d * d * d / 6
L = (
(5 - 2 * c1 + 28 * t1 - 3 * c1 * c1 + 8 * e2lin + 24 * t1 * t1)
* (d ** 5)
* (d**5)
/ 120
)

Expand Down Expand Up @@ -3447,8 +3465,8 @@ def calculateEarthRadius(self, lat, datum):
# Calculate the Earth Radius in meters
eRadius = np.sqrt(
(
(np.cos(lat) * (semiMajorAxis ** 2)) ** 2
+ (np.sin(lat) * (semiMinorAxis ** 2)) ** 2
(np.cos(lat) * (semiMajorAxis**2)) ** 2
+ (np.sin(lat) * (semiMinorAxis**2)) ** 2
)
/ ((np.cos(lat) * semiMajorAxis) ** 2 + (np.sin(lat) * semiMinorAxis) ** 2)
)
Expand Down
Loading