snapshot...

modules/util/resample.py

+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
+