diff --git a/modules/util/resample.py b/modules/util/resample.py index 93d601b943f4a178f00dcc27fd31995cd76cf901..86146634e51c825b38f02eafbc2365b534c8be79 100644 --- a/modules/util/resample.py +++ b/modules/util/resample.py @@ -1,9 +1,11 @@ +import h5py import numpy as np from scipy.interpolate import RegularGridInterpolator, griddata from cartopy.crs import LambertAzimuthalEqualArea import cartopy.crs as ccrs from netCDF4 import Dataset import xarray as xr +from util.util import get_grid_values_all # resample methods: linear = 'linear' @@ -59,7 +61,7 @@ def resample(scalar_field, y_d, x_d, y_t, x_t, method='linear'): return fld_repro -def reproject(fld_2d, lat_2d, lon_2d, proj, target_grid=None, grid_spacing=15000, method=linear): +def reproject(fld_2d, lat_2d, lon_2d, proj, target_grid=None, grid_spacing=15000.0, method=linear): """ :param fld_2d: the 2D scalar field to reproject :param lat_2d: 2D latitude of the scalar field domain @@ -109,3 +111,30 @@ def bisect_great_circle(lon_a, lat_a, lon_b, lat_b): return lon_c, lat_c + +def acspo_sst(filename): + h5f = h5py.File(filename, 'r') + sst = get_grid_values_all(h5f, 'sea_surface_temperature') + sst = sst[0, :, :] + print(sst.shape) + lons = get_grid_values_all(h5f, 'lon') + lats = get_grid_values_all(h5f, 'lat') + print(lons.shape, lats.shape) + + sst = sst[::2, ::2] + lons = lons[::2, ::2] + lats = lats[::2, ::2] + + ylen, xlen = lons.shape + print(ylen, xlen) + + cen_lon = lons[ylen // 2, xlen // 2] + cen_lat = lats[ylen // 2, xlen // 2] + print(cen_lat, cen_lon) + + proj = get_projection('LambertAzimuthalEqualArea', cen_lat, cen_lon) + + fld_reproj, (y_map_2d, x_map_2d) = reproject(sst, lats, lons, proj, grid_spacing=2000.0) + + return fld_reproj, proj, y_map_2d, x_map_2d +