Commit 9b5c8627 authored by Ludvik Brodl's avatar Ludvik Brodl
Browse files

Adds the scripts and setup files

parent 2f49ecf2
# mesan-extractor
Extracts 2018 mesan-data to netcdf4 format.
These scripts were created to create get mesan data into gadget.
\ No newline at end of file
These scripts were created to create get mesan data into gadget.
```
conda create -c conda-forge --yes --name cartopy cartopy python=3.7
conda active cartopy
pip install requirements.txt
```
## Plot mesan data with different projections
```
python mesan_plot.py
```
## Extract mesan data to netcdf4 files. Output is cwd.
```
python mesan2netcdf.py
```
\ No newline at end of file
import datetime
import glob
import itertools
import os
import numpy as np
import netCDF4
import pupygrib
import pyproj
import scipy.interpolate as spint
import scipy.spatial.qhull as qhull
"""
Reprojects and interpolates MESAN 2018 on nsc.bi.liu.se to
SWEREF99TM. Size of output grid is specified by:
OUTGRID_NX
OUTGRID_NY
X_LL, X_UR
Y_LL, Y_UR
Output grid is in regular(reguljära) coordinates in SWEREF99TM. (Square grids)
"""
# Not used, only for documentation:
# ParNr: Name Unit TypeOfLevel height
MESAN_GRIB_TABLE = {
1: ("Pressure", "Pa", "heightAboveSea", 0),
11: ("Temperature", "K", "heightAboveGround", 2),
12: ("Wet bulb temperature", "K", "heightAboveGround", 2),
15: ("Maximum temperature", "K", "heightAboveGround", 2),
16: ("Minimum temperature", "K", "heightAboveGround", 2),
20: ("Visibility", "m", "heightAboveGround", 2),
32: ("Wind gust", "m/s", "heightAboveGround", 10),
33: ("U-component of wind", "m/s", "heightAboveGround", 10),
34: ("V-component of wind", "m/s", "heightAboveGround", 10),
52: ("Relative humidity", "fraction", "heightAboveGround", 2),
71: ("Total cloud cover", "fraction", "heightAboveGround", 0),
73: ("Low cloud cover", "fraction", "heightAboveGround", 0),
74: ("Medium cloud cover", "fraction", "heightAboveGround", 0),
75: ("High cloud cover", "fraction", "heightAboveGround", 0),
77: ("Fraction of significant clouds", "fraction", "heightAboveGround", 0),
78: ("cloud base of significant clouds", "m", "heightAboveSea", 0),
79: ("Cloud Top of significant clouds", "m", "heightAboveGround", 0),
144: ("Frozen part of total precipitation", "fraction", "heightAboveGround", 0),
145: ("Type of precipitation", "code", "heightAboveGround", 0),
146: ("Sort of precipitation", "code", "heightAboveGround", 0),
162: ("12 hour precipitation", "mm", "heightAboveGround", 0),
164: ("24 hour precipitation", "mm", "heightAboveGround", 0),
165: ("1 hour precipitation", "mm", "heightAboveGround", 0),
167: ("3 hour precipitation", "mm", "heightAboveGround", 0),
172: ("12 hour snow", "cm", "heightAboveGround", 0),
174: ("24 hour snow", "cm", "heightAboveGround", 0),
175: ("1 hour snow", "cm", "heightAboveGround", 0),
177: ("3 hour snow", "cm", "heightAboveGround", 0),
}
# All keys in output_netcdf_table will be read from input GRIBs and their values
# specify the output netcdf's parameter.
# 33 and 34 are handled a bit differently since they must first be converted
# from cartesian to polar.
output_netcdf_table = {
1: {"parameter": "P", "long_name": "Pressure", "units": "Pa"},
11: {"parameter": "T2M", "long_name": "Temperature", "units": "K"},
33: {"parameter": "WSPD", "long_name": "Wind speed", "units": "m/s"},
# 33 on input GRIB specifies U-component of wind in m/s
34: {
"parameter": "WDIR",
"long_name": "Wind direction, north: 0 degrees, increasing clockwise",
"units": "Degrees",
}, # 34 on input GRIB specifies V-component of wind in m/s
52: {"parameter": "RH", "long_name": "Relative humidity", "units": "fraction"},
71: {"parameter": "TCC", "long_name": "Total cloud cover", "units": "fraction"},
165: {"parameter": "1H_PREC", "long_name": "1 hour precipitation", "units": "mm"},
175: {"parameter": "1H_SNOW", "long_name": "1 hour snow", "units": "cm"},
}
ROTATED_CRS = pyproj.Proj(
"+proj=ob_tran +o_proj=eqc +lon_0=15 +o_lat_p=30 +a=57.29577953604224884 +b=57.29577953604224884"
)
SWEREF_CRS = pyproj.Proj(init="EPSG:3006")
WGS84_CRS = pyproj.Proj(init="EPSG:4326")
# OUTGRID_X = 740
# OUTGRID_Y = 1570 # 1000m grid
OUTGRID_NX = 74 # 10 km grid
OUTGRID_NY = 157
X_LL, X_UR = (222000, 962000)
Y_LL, Y_UR = (6104000, 7674000)
DX = (X_UR - X_LL) / OUTGRID_NX
DY = (Y_UR - Y_LL) / OUTGRID_NY
OUTGRID_X_COORDS, OUTGRID_Y_COORDS = np.mgrid[X_LL:X_UR:DX, Y_LL:Y_UR:DY]
OUTPUT_GRID_POINTS = np.column_stack(
(list(OUTGRID_X_COORDS.flat), list(OUTGRID_Y_COORDS.flat))
)
OUTGRID_Y_COORDS_AS_WGS84, OUTGRID_X_COORDS_AS_WGS84 = pyproj.transform(
SWEREF_CRS, WGS84_CRS, OUTGRID_Y_COORDS, OUTGRID_X_COORDS
)
OUTGRID_Y_COORDS_AS_WGS84 = OUTGRID_Y_COORDS_AS_WGS84.T
OUTGRID_X_COORDS_AS_WGS84 = OUTGRID_X_COORDS_AS_WGS84.T
GRID_MAPPING_SWEREF99 = {
"false_easting": 500000.0,
"false_northing": 0.0,
"grid_mapping_name": "transverse_mercator",
"latitude_of_projection_origin": 0.0,
"longitude_of_central_meridian": 15.0,
"scale_factor_at_central_meridian": 0.9996,
}
TIME_BOUND_DELTA = 3600 # 1 hour in seconds
# Read, transform and interpolate data of parameters
INPUT_DIR = "/nobackup/smhid14/luftkval/metdata/mesan/2018/"
grib_files = glob.glob(os.path.join(INPUT_DIR, "MESAN_*"))
grib_files.sort()
"""
rotated_latlon_points_in_sweref is extracted from a single
parameter from a single input GRIB file
"""
rotated_latlon_points_in_sweref = None
# Create output netcdf template files without data
template_grib = grib_files[0]
with open(template_grib, "rb") as stream:
for i, msg in enumerate(pupygrib.read(stream), 1):
indicator_of_parameter = msg[1].indicatorOfParameter
if indicator_of_parameter in output_netcdf_table.keys():
print(
f"Setting up output netCDF without data {output_netcdf_table[indicator_of_parameter]}"
)
parameter_name = output_netcdf_table[indicator_of_parameter]["parameter"]
# NetCDF output setup
nc_filename = f"{parameter_name}.nc"
nc_file = netCDF4.Dataset(nc_filename, mode="w", format="NETCDF4")
## Dimensions
nc_file.createDimension("time", None)
nc_file.createDimension("bnds", 2)
nc_file.createDimension("x", OUTGRID_NX)
nc_file.createDimension("y", OUTGRID_NY)
## Create Variables
time_var = nc_file.createVariable("time", "f8", ("time",), zlib=True)
time_bnds_var = nc_file.createVariable(
"time_bnds", "f8", ("time", "bnds"), zlib=True
)
x_var = nc_file.createVariable("x", "f4", ("x"), zlib=True)
y_var = nc_file.createVariable("y", "f4", ("y"), zlib=True)
lon_var = nc_file.createVariable(
"lon", "f8", ("y", "x"), zlib=True, chunksizes=None
)
lat_var = nc_file.createVariable(
"lat", "f8", ("y", "x"), zlib=True, chunksizes=None
)
projection = nc_file.createVariable("EPSG:3006", "S1")
data_var = nc_file.createVariable(
parameter_name, "f4", ("time", "y", "x"), zlib=True, chunksizes=None
)
## Set variables' attributes
### time and time_bnds
time_bnds_var.units = "seconds since 1970-1-1"
time_bnds_var.calendar = "standard"
# times.long_name = "valid time"
time_var.units = "seconds since 1970-1-1"
time_var.bounds = "time_bnds"
time_var.calendar = "standard"
time_var.standard_name = "time"
time_var.axis = "T"
### cf-convention 1.6 grid_mapping
projection.setncatts(GRID_MAPPING_SWEREF99)
### x and y attributes
x_var.standard_name = "projection_x_coordinate"
x_var.units = "m"
x_var.axis = "X"
y_var.standard_name = "projection_y_coordinate"
y_var.units = "m"
y_var.axis = "Y"
### coordinate attributes
lon_var.units = "degrees_east"
lon_var.standard_name = "longitude"
lat_var.units = "degrees_north"
lat_var.standard_name = "latitude"
### Global attributes
nc_file.Conventions = "CF-1.6"
nc_file.institute = (
"Swedish Meteorological and Hydrological Institute (SMHI)"
)
nc_file.history = "This netCDF file created on %s" % datetime.datetime.utcnow().strftime(
"%Y-%m-%d"
)
nc_file.title = "MESAN"
nc_file.source = INPUT_DIR + "MESAN_*"
## Parameters
### get rotated lon lat coordinates
input_lons, input_lats = msg.get_coordinates()
### sweref99 tm x y coordinates
if rotated_latlon_points_in_sweref is None:
x_points, y_points = pyproj.transform(
ROTATED_CRS, SWEREF_CRS, input_lons, input_lats
)
rotated_latlon_points_in_sweref = np.column_stack(
(list(x_points.flat), list(y_points.flat))
)
x_var[:] = OUTGRID_X_COORDS[:, 0]
y_var[:] = OUTGRID_Y_COORDS[0, :]
### transformed lon lat coordinates
lon_var[:, :] = OUTGRID_Y_COORDS_AS_WGS84
lat_var[:, :] = OUTGRID_X_COORDS_AS_WGS84
### data Variable
data_var.grid_mapping = "EPSG:3006"
data_var.coordinates = "lon lat"
data_var.units = output_netcdf_table[indicator_of_parameter]["units"]
data_var.long_name = output_netcdf_table[indicator_of_parameter][
"long_name"
]
# Keep output stream open for later processes
output_netcdf_table[indicator_of_parameter]["netcdf"] = nc_file
output_netcdf_table[indicator_of_parameter]["data_var"] = data_var
output_netcdf_table[indicator_of_parameter]["time_var"] = time_var
output_netcdf_table[indicator_of_parameter]["time_bnds_var"] = time_bnds_var
# Fill output netcdf timesteps and timebounds
START_DT = datetime.datetime(1970, 1, 1)
DATE_FORMAT = "MESAN_%Y%m%d%H%M+000H00M"
grib_files_with_tstep = [] # [(grib_file, tstep, timestamp), ...]
timestamps = []
time_bounds = []
print("\nReading timestamps from GRIB filenames:")
for tstep, grib_file in enumerate(grib_files):
filename = os.path.basename(grib_file)
dt = datetime.datetime.strptime(filename, DATE_FORMAT)
timestamp = int((dt - START_DT).total_seconds())
timestamps.append(timestamp)
time_bounds.append((timestamp, timestamp + TIME_BOUND_DELTA))
grib_files_with_tstep.append((grib_file, tstep))
print(f"Total of {len(timestamps)} timestamps and time bounds.")
print("Adding timestamps and time bounds to output netcdf files:")
for _, output_dict in output_netcdf_table.items():
print(f'{output_dict["netcdf"].filepath()}')
time_var = output_dict["time_var"]
time_bnds_var = output_dict["time_bnds_var"]
time_var[:] = timestamps
time_bnds_var[:] = time_bounds
def uv2wdws(u, v):
"""
Wind component U, V to wind speed and wind direction
"""
a = np.arctan(1.0) * 4.0 / 180.0
wd = 180 + np.arctan2(u, v) / a
ws = np.hypot(u, v)
return ws, wd
def interp_weights(xy, uv, d=2):
"""
Return interpolation weights between 2 irregular grids.
xy: source grid
uv: destination grid
d: Dimension of the grid to be interpolated.
return: interpolation weight matrix
See:
https://stackoverflow.com/questions/20915502/speedup-scipy-griddata-for-multiple-interpolations-between-two-irregular-grids
"""
tri = qhull.Delaunay(xy)
simplex = tri.find_simplex(uv)
vertices = np.take(tri.simplices, simplex, axis=0)
temp = np.take(tri.transform, simplex, axis=0)
delta = uv - temp[:, d]
bary = np.einsum("njk,nk->nj", temp[:, :d, :], delta)
return vertices, np.hstack((bary, 1 - bary.sum(axis=1, keepdims=True)))
def interpolate(values, vtx, wts):
return np.einsum("nj,nj->n", np.take(values, vtx), wts)
vtx, wts = interp_weights(rotated_latlon_points_in_sweref, OUTPUT_GRID_POINTS)
# Fill netcdf files with the parameter data
# 1) Read grid
# 2) Transform grid to sweref points (only done once, see above)
# 3) Interpolate to the output grid using interpolation weights
# from above (linear interpolation).
# 4) Write to netCDF files, 1 timestamp at a time.
START_DT = datetime.datetime(1970, 1, 1)
DATE_FORMAT = "MESAN_%Y%m%d%H%M+000H00M"
for grib_file, tstep in grib_files_with_tstep:
with open(grib_file, "rb") as stream:
print(f"Processing timestep {tstep + 1} / {len(timestamps)} ", end="")
u_wind_component = None
v_wind_component = None
for i, msg in enumerate(pupygrib.read(stream), 1):
indicator_of_parameter = msg[1].indicatorOfParameter
if indicator_of_parameter in output_netcdf_table.keys():
values = msg.get_values()
interpolated_values = interpolate(values, vtx, wts)
reshaped_values = np.reshape(
interpolated_values, (OUTGRID_NX, OUTGRID_NY)
)
if indicator_of_parameter == 33: # U-component
u_wind_component = reshaped_values
elif indicator_of_parameter == 34: # V-component
v_wind_component = reshaped_values
else:
parameter = output_netcdf_table[indicator_of_parameter]["parameter"]
print(parameter, end=" ")
data_var = output_netcdf_table[indicator_of_parameter]["data_var"]
data_var[tstep, :, :] = reshaped_values.T
if (u_wind_component is not None) and (v_wind_component is not None):
wind_speed, wind_dir = uv2wdws(u_wind_component, v_wind_component)
print(output_netcdf_table[33]["parameter"], end=" ")
wspd_data_var = output_netcdf_table[33]["data_var"]
wspd_data_var[tstep, :, :] = wind_speed.T
print(output_netcdf_table[34]["parameter"], end=" ")
wdir_data_var = output_netcdf_table[34]["data_var"]
wdir_data_var[tstep, :, :] = wind_dir.T
# Reset wind component python variables to signal
# complete for this timestamp.
u_wind_component = None
v_wind_component = None
print("")
import matplotlib.pyplot as plt
import pupygrib
import pyproj
import numpy as np
import cartopy.crs as ccrs
from scipy import interpolate
# ParNr: Name Unit TypeOfLevel height
table = {
1: ("Pressure", "Pa", "heightAboveSea", 0),
11: ("Temperature", "K", "heightAboveGround", 2),
12: ("Wet bulb temperature", "K", "heightAboveGround", 2),
15: ("Maximum temperature", "K", "heightAboveGround", 2),
16: ("Minimum temperature", "K", "heightAboveGround", 2),
20: ("Visibility", "m", "heightAboveGround", 2),
32: ("Wind gust", "m/s", "heightAboveGround", 10),
33: ("U-component of wind", "m/s", "heightAboveGround", 10),
34: ("V-component of wind", "m/s", "heightAboveGround", 10),
52: ("Relative humidity", "fraction", "heightAboveGround", 2),
71: ("Total cloud cover", "fraction", "heightAboveGround", 0),
73: ("Low cloud cover", "fraction", "heightAboveGround", 0),
74: ("Medium cloud cover", "fraction", "heightAboveGround", 0),
75: ("High cloud cover", "fraction", "heightAboveGround", 0),
77: ("Fraction of significant clouds", "fraction", "heightAboveGround", 0),
78: ("cloud base of significant clouds", "m", "heightAboveSea", 0),
79: ("Cloud Top of significant clouds", "m", "heightAboveGround", 0),
144: ("Frozen part of total precipitation", "fraction", "heightAboveGround", 0),
145: ("Type of precipitation", "code", "heightAboveGround", 0),
146: ("Sort of precipitation", "code", "heightAboveGround", 0),
162: ("12 hour precipitation", "mm", "heightAboveGround", 0),
164: ("24 hour precipitation", "mm", "heightAboveGround", 0),
165: ("1 hour precipitation", "mm", "heightAboveGround", 0),
167: ("3 hour precipitation", "mm", "heightAboveGround", 0),
172: ("12 hour snow", "cm", "heightAboveGround", 0),
174: ("24 hour snow", "cm", "heightAboveGround", 0),
175: ("1 hour snow", "cm", "heightAboveGround", 0),
177: ("3 hour snow", "cm", "heightAboveGround", 0),
}
with open(
"/nobackup/smhid14/luftkval/metdata/mesan/2018/MESAN_201801010000+000H00M", "rb"
) as stream:
for i, msg in enumerate(pupygrib.read(stream), 1):
lons, lats = msg.get_coordinates()
values = msg.get_values()
print("Message {}: {:.3f} {}".format(i, values.mean(), lons.shape))
lennarts_crs = pyproj.Proj(
"+proj=ob_tran +o_proj=eqc +o_lat_p=30 +lon_0=-10 +a=57.29577953604224884 +b=57.29577953604224884"
)
rotated_crs = pyproj.Proj(
"+proj=ob_tran +o_proj=eqc +lon_0=15 +o_lat_p=30 +a=57.29577953604224884 +b=57.29577953604224884"
)
sweref_crs = pyproj.Proj(init="EPSG:3006")
wgs84_crs = pyproj.Proj(init="EPSG:4326")
print(f"ob_tran: {rotated_crs}")
print(f"wgs84: {wgs84_crs}")
print(f"sweref: {sweref_crs}")
# set up a map
plt.figure(figsize=(6, 7))
rotated_pole = ccrs.RotatedPole(pole_latitude=30, pole_longitude=-165)
# 1st column: Plot GRIB data directly. Cartopy on the fly transformation
ax = plt.subplot(241, projection=ccrs.PlateCarree())
ax.contourf(lons, lats, values, transform=rotated_pole)
ax.set_global()
ax.coastlines()
ax = plt.subplot(245, projection=rotated_pole)
ax.contourf(lons, lats, values, transform=rotated_pole)
ax.set_global()
ax.coastlines()
# 2nd column: Change CRS from rotated latlong to WGS84 datum
# Plot with WGS84 coordinates
lons_transed, lats_transed = pyproj.transform(
rotated_crs, wgs84_crs, lons, lats
)
ax = plt.subplot(242, projection=ccrs.PlateCarree())
ax.contourf(lons_transed, lats_transed, values, transform=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
ax = plt.subplot(246, projection=rotated_pole)
ax.contourf(lons_transed, lats_transed, values, transform=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
# 3rd column: Change CRS to WGS84, Interpolate to a new datum grid (grid_x, grid_y)
# Plot with WGS84 transformed coordinates
grid_x, grid_y = np.mgrid[5:30:100j, 50:70:100j]
lons_transed, lats_transed = pyproj.transform(
rotated_crs, wgs84_crs, lons, lats
)
points = np.column_stack((list(lons_transed.flat), list(lats_transed.flat)))
f = interpolate.griddata(points, list(values.flat), (grid_x, grid_y))
ax = plt.subplot(243, projection=ccrs.PlateCarree())
ax.contourf(grid_x, grid_y, f, transform=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
ax = plt.subplot(247, projection=rotated_pole)
ax.contourf(grid_x, grid_y, f, transform=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
# 4th column: Change CRS to sweref99 and interpolate to geof grid.
# Change CRS to WGS84 and plot in WGS84 transformed coordinates
grid_x, grid_y = np.mgrid[222000:962000:1000, 6104000:7674000:1000]
x, y = pyproj.transform(rotated_crs, sweref_crs, lons, lats)
points = np.column_stack((list(x.flat), list(y.flat)))
f = interpolate.griddata(points, list(values.flat), (grid_x, grid_y))
x_inverse, y_inverse = pyproj.transform(sweref_crs, wgs84_crs, grid_x, grid_y)
ax = plt.subplot(244, projection=ccrs.PlateCarree())
ax.contourf(x_inverse, y_inverse, f, transform=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
ax = plt.subplot(248, projection=rotated_pole)
ax.contourf(x_inverse, y_inverse, f, transform=ccrs.PlateCarree())
ax.set_global()
ax.coastlines()
plt.show()
appdirs==1.4.3
asn1crypto==0.24.0
atomicwrites==1.3.0
attrs==19.1.0
black==19.3b0
bleach==3.1.0
Cartopy==0.17.0
certifi==2019.3.9
cffi==1.12.3
cftime==1.0.3.4
chardet==3.0.4
check-manifest==0.39
Click==7.0
coverage==4.5.3
cryptography==2.7
cycler==0.10.0
docutils==0.14
entrypoints==0.3
filelock==3.0.12
flake8==3.7.7
idna==2.8
importlib-metadata==0.18
jedi==0.13.3
kiwisolver==1.1.0
lxml==4.3.4
matplotlib==3.1.0
mccabe==0.6.1
more-itertools==7.0.0
netCDF4==1.5.1.2
numpy==1.16.4
olefile==0.46
OWSLib==0.17.1
packaging==19.0
parso==0.4.0
Pillow==6.0.0
pkginfo==1.5.0.1
pluggy==0.12.0
pudb==2019.1
-e git+https://notabug.org/mjakob/pupygrib.git@56dd2b33e4db77bde5ebf0f6f0f4d194c892e101#egg=pupygrib
py==1.8.0
pycodestyle==2.5.0
pycparser==2.19
pyepsg==0.4.0
pyflakes==2.1.1
Pygments==2.4.2
pykdtree==1.3.1
pyOpenSSL==19.0.0
pyparsing==2.4.0
pyproj==1.9.6
pyshp==2.1.0
PySocks==1.7.0
pytest==4.6.3
pytest-cov==2.7.1
python-dateutil==2.8.0
pytz==2019.1
readme-renderer==24.0
requests==2.22.0
requests-toolbelt==0.9.1
rope==0.14.0
scipy==1.3.0
Shapely==1.6.4.post2
six==1.12.0
toml==0.10.0
tornado==6.0.2
tox==3.12.1
tqdm==4.32.1
twine==1.13.0
urllib3==1.24.3
urwid==2.0.1
virtualenv==16.6.0
wcwidth==0.1.7
webencodings==0.5.1
zipp==0.5.1
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment