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:
@@ -72,7 +70,7 @@ def reproject(fld_2d, lat_2d, lon_2d, proj, target_grid=None, grid_spacing=15000
            automatically. The target grid is always returned.
     :param grid_spacing: distance between the target grid points (in meters)
     :param method: resampling method: 'linear', 'nearest', 'cubic'
-    :return: reprojected 2D scalar field, the target grid (will be 2D if rotate=True)
+    :return: reprojected 2D scalar field, the target grid
     """
 
     data_xy = proj.transform_points(ccrs.PlateCarree(), lon_2d, lat_2d)[..., :2]
@@ -112,29 +110,37 @@ def bisect_great_circle(lon_a, lat_a, lon_b, lat_b):
     return lon_c, lat_c
 
 
-def acspo_sst(filename):
+def acspo_sst(filename, stride=2, grid_spacing=2000):
     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]
+    # operational
+    # sst = get_grid_values_all(h5f, 'sea_surface_temperature')
+    # sst = sst[0, :, :]
+    # lons = get_grid_values_all(h5f, 'lon')
+    # lats = get_grid_values_all(h5f, 'lat')
+
+    # classic
+    sst = get_grid_values_all(h5f, 'sst_regression')
+    sst = sst[:, :]
+    lons = get_grid_values_all(h5f, 'longitude')
+    lats = get_grid_values_all(h5f, 'latitude')
+
+    print('data dims: ', lons.shape, lats.shape)
+
+    sst = sst[::stride, ::stride]
+    lons = lons[::stride, ::stride]
+    lats = lats[::stride, ::stride]
 
     ylen, xlen = lons.shape
-    print(ylen, xlen)
+    print('final dims: ', ylen, xlen)
 
-    cen_lon = lons[ylen // 2, xlen // 2]
-    cen_lat = lats[ylen // 2, xlen // 2]
-    print(cen_lat, cen_lon)
+    cen_lon = lons[ylen // stride, xlen // stride]
+    cen_lat = lats[ylen // stride, xlen // stride]
+    print('center latitude/longitude: ', 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)
+    fld_reproj, (y_map_2d, x_map_2d) = reproject(sst, lats, lons, proj, grid_spacing=grid_spacing)
 
     return fld_reproj, proj, y_map_2d, x_map_2d