Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Era5 fix. Typo, and need to combine years into single file again #170

Merged
merged 3 commits into from
Apr 30, 2024
Merged

Era5 fix. Typo, and need to combine years into single file again #170

Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
53dcec0
fix renaming typo and combine all outputs into single year rather tha…
ashjbarnes Apr 30, 2024
aca3f48
fix typo and combine all years of era5 into one file
ashjbarnes Apr 30, 2024
171ef0a
black
ashjbarnes Apr 30, 2024
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
145 changes: 73 additions & 72 deletions regional_mom6/regional_mom6.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1572,85 +1572,86 @@ def setup_era5(self, era5_path):
years = [
i for i in range(self.date_range[0].year, self.date_range[1].year + 1)
]
# Loop through each year and read the corresponding files
for year in years:
ds = xr.open_mfdataset(
f"{era5_path}/{fname}/{year}/{fname}*",
decode_times=False,
chunks={"longitude": 100, "latitude": 100},
)
# construct a list of all paths for all years to use for open_mfdataset
paths_per_year = [Path(f"{era5_path}/{fname}/{year}/") for year in years]
all_files = []
for path in paths_per_year:
# Use glob to find all files that match the pattern
files = list(path.glob(f"{fname}*.nc"))
# Add the files to the all_files list
all_files.extend(files)

ds = xr.open_mfdataset(
all_files,
decode_times=False,
chunks={"longitude": 100, "latitude": 100},
)

## Cut out this variable to our domain size
rawdata[fname] = longitude_slicer(
ds,
self.longitude_extent,
"longitude",
).sel(
latitude=slice(
self.longitude_extent[1], self.longitude_extent[0]
) ## This is because ERA5 has latitude in decreasing order (??)
)
## Cut out this variable to our domain size
rawdata[fname] = longitude_slicer(
ds,
self.longitude_extent,
"longitude",
).sel(
latitude=slice(
self.latitude_extent[1], self.latitude_extent[0]
) ## This is because ERA5 has latitude in decreasing order (??)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

What is this question here? Can’t we find out if that’s true or not and remove the “(??)”?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not really sure tbh - might have been a reason for this but I never wrote it down. Don't have too much time today to troubleshoot but was just trying to get this patched for Anupiya (new student of Callum) who couldn't get the example working

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Need to make the most of Gadi being up before it takes a 2 day holiday!!

)

## Now fix up the latitude and time dimensions
## Now fix up the latitude and time dimensions

rawdata[fname] = (
rawdata[fname]
.isel(latitude=slice(None, None, -1)) ## Flip latitude
.assign_coords(
time=np.arange(
0, rawdata[fname].time.shape[0], dtype=float
) ## Set the zero date of forcing to start of run
)
rawdata[fname] = (
rawdata[fname]
.isel(latitude=slice(None, None, -1)) ## Flip latitude
.assign_coords(
time=np.arange(
0, rawdata[fname].time.shape[0], dtype=float
) ## Set the zero date of forcing to start of run
)
)

rawdata[fname].time.attrs = {
"calendar": "julian",
"units": f"hours since {self.date_range[0].strftime('%Y-%m-%d %H:%M:%S')}",
} ## Fix up calendar to match

if fname == "2d":
## Calculate specific humidity from dewpoint temperature
dewpoint = 8.07131 - 1730.63 / (
233.426 + rawdata["2d"]["d2m"] - 273.15
)
humidity = (
(0.622 / rawdata["sp"]["sp"]) * (10**dewpoint) * 101325 / 760
)
q = xr.Dataset(data_vars={"q": humidity})

q.q.attrs = {"long_name": "Specific Humidity", "units": "kg/kg"}
q.to_netcdf(
f"{self.mom_input_dir}/forcing/q_ERA5.nc",
unlimited_dims="time",
encoding={"q": {"dtype": "double"}},
)
elif fname == "crr":
## Calculate total rain rate from convective and total
trr = xr.Dataset(
data_vars={
"trr": rawdata["crr"]["crr"] + rawdata["lsrr"]["lsrr"]
}
)
rawdata[fname].time.attrs = {
"calendar": "julian",
"units": f"hours since {self.date_range[0].strftime('%Y-%m-%d %H:%M:%S')}",
} ## Fix up calendar to match

if fname == "2d":
## Calculate specific humidity from dewpoint temperature
dewpoint = 8.07131 - 1730.63 / (233.426 + rawdata["2d"]["d2m"] - 273.15)
humidity = (0.622 / rawdata["sp"]["sp"]) * (10**dewpoint) * 101325 / 760
q = xr.Dataset(data_vars={"q": humidity})

q.q.attrs = {"long_name": "Specific Humidity", "units": "kg/kg"}
q.to_netcdf(
f"{self.mom_input_dir}/forcing/q_ERA5.nc",
unlimited_dims="time",
encoding={"q": {"dtype": "double"}},
)
elif fname == "crr":
## Calculate total rain rate from convective and total
trr = xr.Dataset(
data_vars={"trr": rawdata["crr"]["crr"] + rawdata["lsrr"]["lsrr"]}
)

trr.trr.attrs = {
"long_name": "Total Rain Rate",
"units": "kg m**-2 s**-1",
}
trr.to_netcdf(
f"{self.mom_input_dir}/forcing/trr_ERA5-{year}.nc",
unlimited_dims="time",
encoding={"trr": {"dtype": "double"}},
)
trr.trr.attrs = {
"long_name": "Total Rain Rate",
"units": "kg m**-2 s**-1",
}
trr.to_netcdf(
f"{self.mom_input_dir}/forcing/trr_ERA5.nc",
unlimited_dims="time",
encoding={"trr": {"dtype": "double"}},
)

elif fname == "lsrr":
## This is handled by crr as both are added together to calculate total rain rate.
pass
else:
rawdata[fname].to_netcdf(
f"{self.mom_input_dir}/forcing/{fname}_ERA5-{year}.nc",
unlimited_dims="time",
encoding={vname: {"dtype": "double"}},
)
elif fname == "lsrr":
## This is handled by crr as both are added together to calculate total rain rate.
pass
else:
rawdata[fname].to_netcdf(
f"{self.mom_input_dir}/forcing/{fname}_ERA5.nc",
unlimited_dims="time",
encoding={vname: {"dtype": "double"}},
)


class segment:
Expand Down
Loading