Einfache, skriptfähige Möglichkeit, unstrukturierte THREDDS-Daten unterprobiern zu können?

Question

Ich versuche, eine Teilmenge von Daten aus einem Dreiecksnetzmodell zu erhalten, das von THREDDS bedient wird. Ich möchte in der Lage sein, eine LAT/LON Bounding Box anzugeben und nur die Daten aus dieser Box zu bekommen. Die Daten-URL ist:

http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_MET_FORECAST.nc

Mit gerasterten Daten es recht einfach ist, die Daten von einem Server THREDDS die Teilmenge. Weiß jemand, was der beste Weg ist, eine Subdomain eines Dreiecksnetzes zu erhalten, das von THREDDS bedient wird?

Für Rasterdaten verwende ich Ferret als meinen OPeNDAP-Client und bin in der Lage, den Download-Prozess zu skripten. Ich möchte hier etwas Ähnliches machen, obwohl ich Matlab, Python oder andere Tools verwenden könnte.

Danke,

Steve

Answer 1

Opendap und NetCDF erlauben keine Extraktion eine unregelmäßige Indizierung. Sie können nur einen Start, Stopp und Schritt anfordern.

Und weil dies ein Dreiecksgitter ist, gibt es keine Garantie, dass Knoten der Dreiecke in der gleichen Region ähnliche Indizes haben. Wenn Sie also nur diese Knoten innerhalb einer Bounding Box erhalten möchten, müssen Sie diese einzeln nacheinander anfordern. Und das ist langsam. In vielen Fällen ist es daher schneller, die minimalen und maximalen Indizes zu bestimmen und den gesamten Chunk in einem Stück anzufordern und dann die Indizes nach Bedarf zu extrahieren.

Hier ist ein Beispielvergleich der beiden Ansätze in Python. In diesem Beispiel wird die Teilmenge zu extrahieren, die alle Indizes umfasst etwa 10-mal schneller als über jeden Index Looping, Zeitreihen zu extrahieren:

import netCDF4 
import time 
import numpy as np 

url='http://www.smast.umassd.edu:8080/thredds/dodsC/FVCOM/NECOFS/Forecasts/NECOFS_GOM3_FORECAST.nc' 
nc = netCDF4.Dataset(url) 
ncv = nc.variables 
lon = ncv['lon'][:] 
lat = ncv['lat'][:] 
tim = ncv['time'][:] 

# find indices inside box 
box = [-71.4,41,-70.2,41.5] 
ii = (lon>=box[0])&(lon<=box[2])&(lat>=box[1])&(lat<=box[3]) 
# jj will have just indices from inside the box: 
jj = np.where(ii)[0] 

Wenn wir eine Schleife über jeden Index, ist es langsam:

# loop over indices, extracting time series 
time0=time.time() 
zi = np.zeros((len(tim),len(jj))) 
k=0 
for j in jj: 
    zi[:,k]=ncv['zeta'][:,j] 
    k +=1 
print('elapsed time: %d seconds' % (time.time()-time0)) 

elapsed time: 56 seconds 

Aber

wenn wir Schleife über den Bereich zu jedem Zeitschritt, es ist viel schneller:

time0=time.time() 
zi2 = np.zeros((len(tim),len(jj))) 
jmin=jj.min() 
jmax=jj.max() 

for i in range(len(tim)): 
    ztmp = ncv['zeta'][i,jmin:jmax+1] 
    zi2[i,:] = ztmp[jj-jmin] 
print('elapsed time: %d seconds' % (time.time()-time0)) 

elapsed time: 6 seconds 

natürlich können Sie die Ergebnisse variieren in Abhängigkeit von der Größe des unstrukturierten Gitter, die Nähe des p In Ihrer Teilmenge finden Sie die Anzahl der Punkte, die Sie extrahieren, usw. Aber das gibt Ihnen hoffentlich eine Vorstellung von den damit verbundenen Problemen.

Answer 2

Richs Antwort bietet die beste Möglichkeit, dies zu tun, wenn Sie sehr enge Leistungsbeschränkungen haben. Aber, wie er sagt, können die Ergebnisse nach Mesh- und Subset-Region variieren.

Wenn Sie daran interessiert ist nur in 1 Zeitschritt zu einem Zeitpunkt (kein Wortspiel beabsichtigt) sind, gibt es Grenzkosten und viel einfacher Code, nur den gesamten Raumbereich greift aus dem dap-Server, wie sie ist:

time0 = time.time() 
zi3 = ncv['zeta'][0, :] 
zi3 = zi3[jj] 
print('elapsed time: %d seconds' % (time.time() - time0)) 

elapsed time: 0 seconds 

Bei der Profilierung (selbst wenn die Arrays vorab zugewiesen werden) ist die grobe min/max-Teilmenge, die Rich verwendet, etwa doppelt so schnell, wenn man dieses Mesh und die spezifische räumliche Teilmenge verwendet. Das NECOFS-Netz ist in der Welt der unstrukturierten Netze relativ klein. Sie könnten zum Beispiel Probleme mit ADCIRC-Netzen im atlantischen Maßstab bekommen.

POSTSCRIPT: Sie greifen das gesamte Gitter, wenn Sie lat und lon sowieso, um Ihre Subsetting-Indizes zu bestimmen.