Przeprowadź pętlę przez pliki netcdf i wykonaj obliczenia - Python lub R

Question

To jest mój pierwszy raz, kiedy używam netCDF i staram się owijać moją głowę.

Mam wiele wersji 3 plików netcdf (NOAA NARR air.2m średnie dzienne przez cały rok). Każdy plik obejmuje rok od 1979 do 2012 roku. Są to siatki 349 x 277 o rozdzielczości około 32 km. Dane zostały pobrane z here.

Wymiar to czas (godziny od 1 stycznia 1800), a moją zmienną będącą przedmiotem zainteresowania jest powietrze. Trzeba obliczyć skumulowane dni z temperaturą < 0. Na przykład

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 

muszę przechowywać te dane w nowym pliku netcdf. Jestem zaznajomiony z Pythonem i nieco z R. Jak najlepiej przechodzić przez każdy dzień, sprawdzać wartość poprzednich dni i na tej podstawie wyprowadzać wartość do nowego pliku netcdf z tym samym wymiarem i zmienną ... lub może po prostu dodaj inną zmienną do oryginalnego pliku netcdf z danymi wyjściowymi, których szukam.

Czy najlepiej pozostawić wszystkie pliki oddzielnie lub połączyć je? Połączyłem je z ncrcat i działało dobrze, ale plik ma 2,3 gb.

Dzięki za wejście.

Mój obecny postęp w Pythonie:

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 

Aby pomóc mi to lepiej zrozumieć, jakiego rodzaju struktury danych mam pracy z? Czy ['air'] klucz w powyższym przykładzie i [1,1,1] są również kluczami? aby uzyskać wartość 298,37473. Jak mogę przejść przez [1,1,1]?

Answer 1

Możesz użyć bardzo ładnej funkcji MFDataset w netCDF4, aby traktować kilka plików jako jeden zagregowany plik, bez potrzeby używania ncrcat. Więc kod będzie wyglądać następująco:

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() 

a oto jeden sposób, aby zapisać plik (nie może być prostsze sposoby):

# 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() 

Gdyby spróbować patrząc na otrzymaną plik w the Unidata NetCDF-Java Tools-UI GUI, wydaje się być w porządku: Zauważ, że tutaj właśnie pobrałem dwa zestawy danych do testowania, więc użyłem

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

jako przykład. Dla wszystkich danych, można użyć

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

lub

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

Answer 2

Oto rozwiązanie R.

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

Wiem, że to dość późno dla tego wątku z 2013 roku, ale po prostu chcę podkreślić, że przyjęte rozwiązanie nie zapewnia rozwiązanie dokładne pytanie postawione.Wydaje się, że pytanie dotyczy długości każdego ciągłego okresu temperatur poniżej zera (uwaga w pytaniu licznika resetuje się, gdy temperatura przekracza zero), co może mieć znaczenie dla zastosowań klimatycznych (np. Dla rolnictwa), podczas gdy przyjęte rozwiązanie daje jedynie sumę liczba dni w roku, w których temperatura jest niższa niż zero. Jeśli tak jest naprawdę to, co mkmitchell chce (to zostało zaakceptowane jako odpowiedź), to można to zrobić z wiersza polecenia w cdo bez konieczności martwienia się o netcdf wejścia/wyjścia:

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

więc zapętlony Skrypt ten być:

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 

Jeśli następnie chce całkowitą liczbę w całym okresie można następnie plików zamiast kot _numdays które są znacznie mniejsze:

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

ale znowu pytanie zobaczyć ms, aby zebrać dni za zdarzenie, co jest inne, ale nie zostało udzielone przez zaakceptowaną odpowiedź.