Schleife durch Netcdf-Dateien und führen Sie Berechnungen - Python oder R

Question

Dies ist mein erstes Mal mit netCDF und ich versuche, meinen Kopf um mich damit zu arbeiten.

Ich habe mehrere Netcdf-Dateien der Version 3 (NOAA NARR air.2m tägliche Mittelwerte für ein ganzes Jahr). Jede Datei erstreckt sich über ein Jahr zwischen 1979 - 2012. Sie sind 349 x 277 Gitter mit einer Auflösung von etwa 32 km. Daten wurden von here heruntergeladen.

Die Dimension ist Zeit (Stunden seit 1/1/1800) und meine Variable von Interesse ist Luft. Ich muss akkumulierte Tage mit einer Temperatur

Day 1 = +4 degrees, accumulated days = 0 
    Day 2 = -1 degrees, accumulated days = 1 
    Day 3 = -2 degrees, accumulated days = 2 
    Day 4 = -4 degrees, accumulated days = 3 
    Day 5 = +2 degrees, accumulated days = 0 
    Day 6 = -3 degrees, accumulated days = 1 

ich diese Daten in einer neuen netcdf Datei speichern müssen < 0. Für

Beispiel berechnen. Ich bin vertraut mit Python und etwas mit R. Was ist der beste Weg, um jeden Tag durchlaufen, überprüfen Sie die vorherigen Tage Wert, und darauf basierend, geben Sie einen Wert in eine neue netcdf-Datei mit genau der gleichen Dimension und Variable ... Oder fügen Sie einfach eine weitere Variable zur ursprünglichen netcdf-Datei mit der Ausgabe hinzu, die ich suche.

Ist es am besten, alle Dateien zu trennen oder zu kombinieren? Ich habe sie mit ncrcat kombiniert und es hat gut funktioniert, aber die Datei ist 2,3 GB.

Danke für die Eingabe.

Mein aktueller Fortschritt in Python:

import numpy 
import netCDF4 
#Change my working DIR 
f = netCDF4.Dataset('air7912.nc', 'r') 
for a in f.variables: 
    print(a) 

#output = 
    lat 
    long 
    x 
    y 
    Lambert_Conformal 
    time 
    time_bnds 
    air 

f.variables['air'][1, 1, 1] 
#Output 
    298.37473 

mir zu helfen, diese besser zu verstehen, welche Art von Datenstruktur bin ich mit zu arbeiten? Ist [Luft] der Schlüssel im obigen Beispiel und [1,1,1] sind auch Schlüssel? um den Wert von 298.37473 zu erhalten. Wie kann ich dann [1,1,1] durchlaufen?

Answer 1

Sie können die sehr nützliche Funktion MFDataset in netCDF4 verwenden, um eine Gruppe von Dateien als eine aggregierte Datei zu behandeln, ohne dass Sie ncrcat verwenden müssen. So würden Sie Code wie folgt aussehen:

from pylab import * 
import netCDF4 

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.19??.nc') 
# print variables 
f.variables.keys() 

atemp = f.variables['air'] 
print atemp 

ntimes, ny, nx = shape(atemp) 
cold_days = zeros((ny,nx),dtype=int) 

for i in xrange(ntimes): 
    cold_days += atemp[i,:,:].data-273.15 < 0 

pcolormesh(cold_days) 
colorbar() 

Und hier ist ein Weg, um die Datei zu schreiben (es könnte einfacher Weise sein):

# create NetCDF file 
nco = netCDF4.Dataset('/usgs/data2/notebook/cold_days.nc','w',clobber=True) 
nco.createDimension('x',nx) 
nco.createDimension('y',ny) 

cold_days_v = nco.createVariable('cold_days', 'i4', ('y', 'x')) 
cold_days_v.units='days' 
cold_days_v.long_name='total number of days below 0 degC' 
cold_days_v.grid_mapping = 'Lambert_Conformal' 

lono = nco.createVariable('lon','f4',('y','x')) 
lato = nco.createVariable('lat','f4',('y','x')) 
xo = nco.createVariable('x','f4',('x')) 
yo = nco.createVariable('y','f4',('y')) 
lco = nco.createVariable('Lambert_Conformal','i4') 

# copy all the variable attributes from original file 
for var in ['lon','lat','x','y','Lambert_Conformal']: 
    for att in f.variables[var].ncattrs(): 
     setattr(nco.variables[var],att,getattr(f.variables[var],att)) 

# copy variable data for lon,lat,x and y 
lono[:]=f.variables['lon'][:] 
lato[:]=f.variables['lat'][:] 
xo[:]=f.variables['x'][:] 
yo[:]=f.variables['y'][:] 

# write the cold_days data 
cold_days_v[:,:]=cold_days 

# copy Global attributes from original file 
for att in f.ncattrs(): 
    setattr(nco,att,getattr(f,att)) 

nco.Conventions='CF-1.6' 
nco.close() 

Wenn ich auf der resultierende versuchen suchen Datei in the Unidata NetCDF-Java Tools-UI GUI, scheint es in Ordnung zu sein: Beachten Sie auch, dass ich hier nur zwei der Datensätze zum Testen heruntergeladen, so dass ich

verwendet

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.19??.nc') 

als ein Beispiel. Für alle die Daten, könnten Sie

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.????.nc') 

oder

f = netCDF4.MFDataset('/usgs/data2/rsignell/models/ncep/narr/air.2m.*.nc') 

Answer 2

Hier eine R Lösung verwenden.

infiles <- list.files("data", pattern = "nc", full.names = TRUE, include.dirs = TRUE) 

outfile <- "data/air.colddays.nc"  

library(raster) 

r <- raster::stack(infiles) 
r <- sum((r - 273.15) < 0) 

plot(r) 

Answer 3

Ich weiß, das ab 2013 für diesen Thread ziemlich spät ist, aber ich möchte nur darauf hinweisen, dass die akzeptierte Lösung nicht die Lösung gestellt erbringt die genaue Frage.Die Frage scheint die Länge jeder kontinuierlichen Periode von Temperaturen unter Null zu sein (beachte die Frage, ob der Zähler zurückgesetzt wird, wenn die Temperatur Null übersteigt), was für Klimaanwendungen (z. B. für die Landwirtschaft) wichtig sein kann, während die akzeptierte Lösung nur die Summe angibt Anzahl der Tage im Jahr, an denen die Temperatur unter Null liegt. Ist dies wirklich das, was mkmitchell will (es als Antwort akzeptiert wurde), dann kann es in in cdo von der Kommandozeile ausgeführt werden, ohne sich um Ein-/Ausgabe netcdf kümmern:

cdo timsum -lec,273.15 in.nc out.nc 

so ein geschlungenen Skript würde sein:

files=`ls *.nc` # pick up all the netcdf files in a directory 
for file in $files ; do 
    # I use 273.15 as from the question seems T is in Kelvin 
    cdo timsum -lec,273.15 $file ${file%???}_numdays.nc 
done 

Wenn Sie dann die Gesamtzahl über den gesamten Zeitraum möchten, können Sie dann die _numdays Katze stattdessen Dateien, die viel kleiner sind:

cdo cat *_numdays.nc total.nc 
cdo timsum total.nc total_below_zero.nc 

Aber auch hier ist die Frage sehen ms wollen akkumulierte Tage haben pro Ereignis, was anders ist, aber nicht durch die angenommene Antwort zur Verfügung gestellt.