Skip to content
Extraits de code Groupes Projets
Valider 2d2eb518 rédigé par Benjamin Richaud's avatar Benjamin Richaud
Parcourir les fichiers

Merge branch 'master' into 'main'

Some cleaning of the code

See merge request !3
parents de8def84 b3bb16e5
Aucune branche associée trouvée
Aucune étiquette associée trouvée
1 requête de fusion!3Some cleaning of the code
......@@ -5,7 +5,6 @@ import matplotlib.path as mpath
import cartopy.crs as ccrs
import cartopy.feature as cfeat
from pathlib import Path
from pyproj import CRS, Transformer
def define_sectors(ref="RH"):
"""
......@@ -86,8 +85,9 @@ def interp_mask2grid(target_grid=None,
Inputs: - target_grid: the target grid on which to interpolate the NSIDC mask.
It can be a xr.Dataset containing lat and lon variables defining the grid,
or a string mentioning the grid:
if weights='ORCA1_nh', use the weights for Northern Hemisphere ORCA1 config;
if 'ORCA025_nh', use the weights for NorthH for ORCA025
- if 'ORCA1_nh', use the mask for Northern Hemisphere ORCA1 config;
- if 'ORCA025_nh', use the mask for Northern Hemisphere for ORCA025;
- any other string: a user-calculated mask already interpolated on the target grid.
#TODO: need to add full ORCA1 and ORCA025, as well as ORCA12 and ORCA12_nh
-maskFile: path and name of file containing the NSIDC mask.
......@@ -101,41 +101,22 @@ def interp_mask2grid(target_grid=None,
ds_mask_interp = xr.open_dataset(Path.cwd().joinpath('../data/NSIDCRegionsMask_ORCA025_nh.nc'))
else:
ds_mask_interp = xr.open_dataset(target_grid)
else:
else: # If the target grid is a dataset or dataarray, interpolate the mask onto it.
# First of all, load the mask
ds_mask = xr.open_dataset(maskFile)
# ------------- If no lat-lon in the ds_mask, calculate and add them -------------
if 'lat' not in list(ds_mask.coords):
# Convert from projected coordinates to lat-lon coordinates
# A first issue is that the mask is provided in North Polar Stereographic projection
# Extract the x and y (projection) coordinates and grid them (from 1D to 2D)
x_mask, y_mask = np.meshgrid(ds_mask.x, ds_mask.y)
# Define the projection and geographic systems according to CRS conventions.
orig_crs = CRS.from_epsg(3413) # North Polar Stereographic from NSIDC
target_crs = CRS.from_epsg(4326) # Geographic coordinate system (i.e. lat-lon)
# Use pyproj.Transformer to convert from NSIDC projection to lat-lon coordinates
transformer = Transformer.from_crs(orig_crs, target_crs)
lat_src, lon_src = transformer.transform(x_mask, y_mask)
# Assign the new coordinates to the dataset
ds_mask = ds_mask.assign_coords({'lon': (['y', 'x'], lon_src), 'lat': (['y', 'x'], lat_src)})
# ------------- Bulk of the function: interpolate the mask to the ds grid -------------
# If 'weights' are given, load them directly. If not, calculate them
# But since xesmf does not exist yet on `cyclone`, need to deal with it.
# try: # If xesmf exists, create or load the regridder and regrid
import xesmf as xe # type: ignore
# Import xesmf. Will throw an error if not available
import xesmf as xe
# Create a regridder.
Regridder = xe.Regridder(ds_mask, target_grid, method='bilinear')
# Now apply this regridder to the mask
ds_mask_interp = Regridder(ds_mask)
# except ImportError: # If the import failed, send an error
# print("ERROR: xesmf could not be loaded, go to `coriolis`.")
# return
# return the interpolated mask.
# Return the interpolated mask.
return ds_mask_interp
def groupby_sectors(ds, ref=None,
target_gd=None,
target_gd=None,
maskFile=Path.cwd().joinpath('../data/NSIDCRegions_N3.125km_v1.1_wLatLon_df1.nc')):
"""
Define sectors and Groupby a dataset or dataArray into those sectors.
......
Ce diff est replié.
0% Chargement en cours ou .
You are about to add 0 people to the discussion. Proceed with caution.
Terminez d'abord l'édition de ce message.
Veuillez vous inscrire ou vous pour commenter