Regración de datos regulares de netcdf

Tengo un archivo netcdf que contiene temperaturas globales de la superficie del mar. Usando matplotlib y Basemap, he logrado hacer un mapa de estos datos, con el siguiente código:

from netCDF4 import Dataset import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.basemap import Basemap filename = '/Users/Nick/Desktop/SST/SST.nc' fh = Dataset(filename, mode='r') lons = fh.variables['LON'][:] lats = fh.variables['LAT'][:] sst = fh.variables['SST'][:].squeeze() fig = plt.figure() m = Basemap(projection='merc', llcrnrlon=80.,llcrnrlat=-25.,urcrnrlon=150.,urcrnrlat=25.,lon_0=115., lat_0=0., resolution='l') lon, lat = np.meshgrid(lons, lats) xi, yi = m(lon, lat) cs = m.pcolormesh(xi,yi,sst, vmin=18, vmax=32) m.drawmapboundary(fill_color='0.3') m.fillcontinents(color='0.3', lake_color='0.3') cbar = m.colorbar(cs, location='bottom', pad="10%", ticks=[18., 20., 22., 24., 26., 28., 30., 32.]) cbar.set_label('January SST (' + u'\u00b0' + 'C)') plt.savefig('SST.png', dpi=300) 

El problema es que los datos son de muy alta resolución (cuadrícula de 9 km), lo que hace que la imagen resultante sea bastante ruidosa. Me gustaría colocar los datos en una cuadrícula de resolución más baja (por ejemplo, 1 grado), pero estoy luchando para averiguar cómo se puede hacer esto. Seguí una solución trabajada para probar y usar la función matplotlib griddata insertando el código a continuación en mi ejemplo anterior, pero resultó en que ‘ValueError: condition debe ser una matriz 1-d’.

 xi, yi = np.meshgrid(lons, lats) X = np.arange(min(x), max(x), 1) Y = np.arange(min(y), max(y), 1) Xi, Yi = np.meshgrid(X, Y) Z = griddata(xi, yi, z, Xi, Yi) 

Soy un principiante relativo de Python y matplotlib, por lo que no estoy seguro de lo que estoy haciendo mal (o de lo que podría ser un mejor enfoque). Cualquier consejo apreciado!

Si vuelve a colocar los datos en una cuadrícula de lat / lon más gruesa utilizando, por ejemplo, interpolación bilineal, esto dará como resultado un campo más suave .

La guía ClimateData de NCAR tiene una buena introducción a la redistribución (en general, no específica de Python).

La implementación más poderosa de las rutinas de regridding disponibles para Python es, por lo que sé, la interfaz de Python (ESMPy) del Framework de Modelado del Sistema Terrestre (ESMF) . Si esto es un poco demasiado complicado para su aplicación, debería considerar

  1. Tutoriales de EarthPy sobre regridding (por ejemplo, usando Pyresample , cKDTree o Basemap ).
  2. Convirtiendo sus datos en un cubo Iris y usando las funciones de redistribución de Iris .

Tal vez empiece por mirar el tutorial de regreding de EarthPy usando el mapa base , ya que ya lo está utilizando.

La forma de hacer esto en tu ejemplo sería

 from mpl_toolkits import basemap from netCDF4 import Dataset filename = '/Users/Nick/Desktop/SST/SST.nc' with Dataset(filename, mode='r') as fh: lons = fh.variables['LON'][:] lats = fh.variables['LAT'][:] sst = fh.variables['SST'][:].squeeze() lons_sub, lats_sub = np.meshgrid(lons[::4], lats[::4]) sst_coarse = basemap.interp(sst, lons, lats, lons_sub, lats_sub, order=1) 

Esto realiza una interpolación bilineal ( order=1 ) en sus datos de SST en una cuadrícula submuestreada (cada cuarto punto). Su ttwig se verá más de grano grueso después. Si no te gusta eso, interpola de nuevo en la cuadrícula original con, por ejemplo,

 sst_smooth = basemap.interp(sst_coarse, lons_sub[0,:], lats_sub[:,0], *np.meshgrid(lons, lats), order=1) 

Normalmente ejecuto mis datos a través de un filtro de Laplace para suavizar. Tal vez podría probar la siguiente función y ver si ayuda con sus datos. La función se puede llamar con o sin una máscara (por ejemplo, máscara de tierra / océano para puntos de datos del océano). Espero que esto ayude. T

 # Laplace filter for 2D field with/without mask # M = 1 on - cells used # M = 0 off - grid cells not used # Default is without masking import numpy as np def laplace_X(F,M): jmax, imax = F.shape # Add strips of land F2 = np.zeros((jmax, imax+2), dtype=F.dtype) F2[:, 1:-1] = F M2 = np.zeros((jmax, imax+2), dtype=M.dtype) M2[:, 1:-1] = M MS = M2[:, 2:] + M2[:, :-2] FS = F2[:, 2:]*M2[:, 2:] + F2[:, :-2]*M2[:, :-2] return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F) def laplace_Y(F,M): jmax, imax = F.shape # Add strips of land F2 = np.zeros((jmax+2, imax), dtype=F.dtype) F2[1:-1, :] = F M2 = np.zeros((jmax+2, imax), dtype=M.dtype) M2[1:-1, :] = M MS = M2[2:, :] + M2[:-2, :] FS = F2[2:, :]*M2[2:, :] + F2[:-2, :]*M2[:-2, :] return np.where(M > 0.5, (1-0.25*MS)*F + 0.25*FS, F) # The mask may cause laplace_X and laplace_Y to not commute # Take average of both directions def laplace_filter(F, M=None): if M == None: M = np.ones_like(F) return 0.5*(laplace_X(laplace_Y(F, M), M) + laplace_Y(laplace_X(F, M), M)) 

Para responder a su pregunta original con respecto a scipy.interpolate.griddata , también:

Observe detenidamente las especificaciones de los parámetros para esa función (por ejemplo, en la documentación de SciPy ) y asegúrese de que sus matrices de entrada tengan las formas correctas . Es posible que tengas que hacer algo como

 import numpy as np points = np.vstack([a.flat for a in np.meshgrid(lons,lats)]).T # (n,D) values = sst.ravel() # (n) 

etc.