open - netcdf in python

Regulando datos regulares de netcdf

Question

Tengo un archivo netcdf que contiene las temperaturas globales de la superficie del mar. Usando matplotlib y Basemap, logré 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 tienen una resolución muy alta (cuadrícula de 9 km) lo que hace que la imagen resultante sea bastante ruidosa. Me gustaría poner los datos en una cuadrícula de resolución más baja (por ejemplo, 1 grado), pero estoy luchando para encontrar la forma de hacerlo. Seguí una solución trabajada para probar y usar la función de cuadrícula de matplotlib insertando el código a continuación en mi ejemplo anterior, pero resultó en ''ValueError: la condición debe ser una matriz de 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 qué mejor enfoque podría tener). Cualquier consejo apreciado!

Answer 1

Por lo general, ejecuto mis datos a través de un filtro de Laplace para suavizar. Tal vez podrías probar la función a continuación y ver si te ayuda con tus datos. La función se puede llamar con o sin una máscara (por ejemplo, máscara terrestre / oceánica para puntos de datos oceánicos). 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))

Answer 2

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

Eche un vistazo de cerca a las especificaciones de 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 deba 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.

Answer 3

Si cambia los datos a una grilla de latitud / longitud más gruesa utilizando, por ejemplo, interpolación bilineal, se obtendrá un campo más uniforme.

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

La implementación más poderosa de las rutinas de actualización disponibles para Python es, que yo sepa, la interfaz Python del Marco de Modelado del Sistema Tierra (ESMF) (ESMPy) . Si esto es demasiado complicado para su aplicación, debe considerar

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

Tal vez comience mirando el tutorial de actualización de EarthPy usando Basemap , ya que lo está usando.

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 la interpolación bilineal ( order=1 ) en sus datos de SST en una cuadrícula submuestreada (cada cuarto punto). Su trama se verá más gruesa después. Si no te gusta, 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)