Lecture de séries chronologiques à partir de netCDF avec python

Question

J'essaie de créer des séries temporelles à partir d'un fichier netCDF (accessible via le serveur Thredds) avec python. Le code que j'utilise semble correct, mais les valeurs de la variable amb lecture sont "masquées". Je suis novice en python et je ne connais pas les formats. Une idée de comment puis-je lire les données?

C'est le code que j'utilise:

import netCDF4 
import pandas as pd 
import datetime as dt 
import matplotlib.pyplot as plt 
from datetime import datetime, timedelta # 

dayFile = datetime.now() - timedelta(days=1) 
dayFile = dayFile.strftime("%Y%m%d") 

url='http://nomads.ncep.noaa.gov:9090/dods/nam/nam%s/nam1hr_00z' %(dayFile) 

# NetCDF4-Python can open OPeNDAP dataset just like a local NetCDF file 
nc = netCDF4.Dataset(url) 
varsInFile = nc.variables.keys() 

lat = nc.variables['lat'][:] 
lon = nc.variables['lon'][:] 
time_var = nc.variables['time'] 
dtime = netCDF4.num2date(time_var[:],time_var.units) 

first = netCDF4.num2date(time_var[0],time_var.units) 
last = netCDF4.num2date(time_var[-1],time_var.units) 
print first.strftime('%Y-%b-%d %H:%M') 
print last.strftime('%Y-%b-%d %H:%M') 

# determine what longitude convention is being used 
print lon.min(),lon.max() 

# Specify desired station time series location 
# note we add 360 because of the lon convention in this dataset 
#lati = 36.605; loni = -121.85899 + 360. # west of Pacific Grove, CA 
lati = 41.4; loni = -100.8 +360.0 # Georges Bank 

# Function to find index to nearest point 
def near(array,value): 
idx=(abs(array-value)).argmin() 
return idx 

# Find nearest point to desired location (no interpolation) 
ix = near(lon, loni) 
iy = near(lat, lati) 
print ix,iy 

# Extract desired times. 

# 1. Select -+some days around the current time: 
start = netCDF4.num2date(time_var[0],time_var.units) 
stop = netCDF4.num2date(time_var[-1],time_var.units) 
time_var = nc.variables['time'] 
datetime = netCDF4.num2date(time_var[:],time_var.units) 

istart = netCDF4.date2index(start,time_var,select='nearest') 
istop = netCDF4.date2index(stop,time_var,select='nearest') 
print istart,istop 

# Get all time records of variable [vname] at indices [iy,ix] 
vname = 'dswrfsfc' 

var = nc.variables[vname] 

hs = var[istart:istop,iy,ix] 

tim = dtime[istart:istop] 

# Create Pandas time series object 
ts = pd.Series(hs,index=tim,name=vname) 

Les données var ne sont pas lues comme je m'y attendais, apparemment parce que les données sont masquées:

>>> hs 
masked_array(data = [-- -- -- ..., -- -- --], 
      mask = [ True True True ..., True True True], 
     fill_value = 9.999e+20) 

Le nom var, et les séries chronologiques sont correct, aussi bien du reste du script. La seule chose qui ne fonctionne pas est la var data récupérée. C'est la série de temps que j'obtiens:

>>> ts 
2016-10-25 00:00:00.000000 NaN 
2016-10-25 01:00:00.000000 NaN 
2016-10-25 02:00:00.000006 NaN 
2016-10-25 03:00:00.000000 NaN 
2016-10-25 04:00:00.000000 NaN 
... ... ... ... ... 
2016-10-26 10:00:00.000000 NaN 
2016-10-26 11:00:00.000006 NaN 
Name: dswrfsfc, dtype: float32 

Toute aide sera appréciée!

Answer 1

Hmm, ce code semble familier. ;-)

Vous obtenez des NaN parce que le modèle NAM auquel vous essayez d'accéder utilise maintenant la longitude dans la plage [-180, 180] au lieu de la plage [0, 360]. Donc, si vous demandez loni = -100.8 au lieu de loni = -100.8 +360.0, je crois que votre code retournera des valeurs non-NaN. Cependant, il est intéressant de noter que la tâche d'extraction de séries temporelles à partir de données multidimensionnelles est maintenant beaucoup plus facile avec xarray, car vous pouvez simplement sélectionner un ensemble de données le plus proche d'un point de latitude, puis tracer n'importe quelle variable. Les données ne sont chargées que lorsque vous en avez besoin, et non lorsque vous extrayez l'objet de l'ensemble de données. Donc, fondamentalement, vous avez besoin maintenant seulement:

import xarray as xr 

ds = xr.open_dataset(url) # NetCDF or OPeNDAP URL 
lati = 41.4; loni = -100.8 # Georges Bank 

# Extract a dataset closest to specified point 
dsloc = ds.sel(lon=loni, lat=lati, method='nearest') 

# select a variable to plot 
dsloc['dswrfsfc'].plot() 

ordinateur portable complet ici: http://nbviewer.jupyter.org/gist/rsignell-usgs/d55b37c6253f27c53ef0731b610b81b4