Rund lat/lon Ernte einer NetCDF-Datei mit Python

Question

Ich bin an einem Projekt arbeiten, die Raw-Binary-Radardaten von der National Weather Service FTP-Site auf einen Server importiert. Mit einem Datenexport-Tool von Weather and Climate Toolkit konvertiere ich die Daten in eine netCDF-Datei. Hier finden Sie das Ergebnis eines „ncdump -h“ Befehl auf der .nc Datei:

netcdf last { 
dimensions: 
     lat = 800 ; 
     lon = 1200 ; 
     time = 1 ; 
variables: 
     double cref(time, lat, lon) ; 
       cref:long_name = "Level-III Composite Reflectivity (16 levels/248 nm)" ; 
       cref:missing_value = -999. ; 
       cref:units = "dBZ" ; 
     double lat(lat) ; 
       lat:units = "degrees_north" ; 
       lat:spacing = "0.010995604400775773" ; 
       lat:datum = "NAD83 - NOAA Standard" ; 
     double lon(lon) ; 
       lon:units = "degrees_east" ; 
       lon:spacing = "0.010983926942902655" ; 
       lon:datum = "NAD83 - NOAA Standard" ; 
     int time(time) ; 
       time:units = "seconds since 1970-1-1" ; 

// global attributes: 
       :title = "Level-III Composite Reflectivity (16 levels/248 nm) 22:23:47 UTC 10/20/2016" ; 
       :Conventions = "CF-1.0" ; 
       :History = "Exported to NetCDF-3 CF-1.0 conventions by the NOAA Weather and Climate Toolkit (version 3.7.9) \n", 
        "Export Date: Thu Oct 20 16:11:07 EDT 2016" ; 
       :geographic_datum_ESRI_PRJ = "GEOGCS[\"GCS_North_American_1983\",DATUM[\"D_North_American_1983\",SPHEROID[\"GRS_1980\",6378137,298.257222101]],PRIMEM[\"Greenwich\",0],UNIT[\"Degree\",0.0174532925199433]]" ; 
       :geographic_datum_OGC_WKT = "GEOGCS[\"NAD83\", DATUM[\"NAD83\", SPHEROID[\"GRS_1980\", 6378137.0, 298.25722210100002],TOWGS84[0,0,0,0,0,0,0]], PRIMEM[\"Greenwich\", 0.0], UNIT[\"degree\",0.017453292519943295], AXIS[\"Longitude\",EAST], AXIS[\"Latitude\",NORTH]]" ; 
} 

Ich möchte den größten Eintrag für die cref Variable finden, die ich mit dem netCDF4 ziemlich leicht tun können, und numpy Bibliotheken in python:

import netCDF4 
import numpy 

netcdf = netCDF4.Dataset("last.nc") 
var = netcdf.variables['cref'] 
print(numpy.nanmax(var)) 
print(numpy.nanmin(var)) 

aber ich hoffe, die netCDF-Dateien, so dass der max und min zu filtern, dass nur in einem gewissen Abstand von einem gegebenen lat/lon zu finden sind. Mit anderen Worten, ich bin die Hoffnung, ein Kreis mit einem bestimmten Radius auf „crop“ um einen bestimmten lat/lon. Ich habe how to crop a square durch einen anderen SO Thread gefunden, kann aber nicht herausfinden, wie ein Kreis funktionieren würde.

Answer 1

würde ich den Abstand zwischen der Mitte und jeder lat/lon Paar (2D-Gitter) berechnen, und verwenden, die eine Maske zu konstruieren, die Sie auf Ihre Daten anwenden können. Sobald maskiert, können Sie wieder einfach numpy Funktionen verwenden, um Statistiken wie max() zu berechnen.

Zum Beispiel mit der haversine() Funktion von https://stackoverflow.com/a/4913653/3581217, auf eine vektorisierte Version geändert, die Sie direkt auf numpy Arrays anwenden können:

import numpy as np 
import matplotlib.pylab as pl 

def haversine(lon1, lat1, lon2, lat2): 
    # convert decimal degrees to radians 
    lon1 = np.deg2rad(lon1) 
    lon2 = np.deg2rad(lon2) 
    lat1 = np.deg2rad(lat1) 
    lat2 = np.deg2rad(lat2) 

    # haversine formula 
    dlon = lon2 - lon1 
    dlat = lat2 - lat1 
    a = np.sin(dlat/2)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2 
    c = 2 * np.arcsin(np.sqrt(a)) 
    r = 6371 
    return c * r 

# Latitude/longitude grid 
lat = np.linspace(50,54,16) 
lon = np.linspace(6,9,12) 

# Center coordinates 
clat = 52 
clon = 7 

max_dist = 100  # max distance in km 

# Calculate distance between center and all other lat/lon pairs 
distance = haversine(lon[:,np.newaxis], lat, clon, clat) 

# Mask distance array where distance > max_dist 
distance_m = np.ma.masked_greater(distance, max_dist) 

# Dummy data 
data = np.random.random(size=[lon.size, lat.size]) 

# Test: set a value outside the max_dist circle to a large value: 
data[0,0] = 10 

# Mask the data array based on the distance mask 
data_m = np.ma.masked_where(distance > max_dist, data) 

pl.figure() 
pl.subplot(221) 
pl.title('distance (km)') 
pl.pcolormesh(lon, lat, np.transpose(distance)) 
pl.colorbar() 

pl.subplot(222) 
pl.title('distance < max_dist (km)') 
pl.pcolormesh(lon, lat, np.transpose(distance_m)) 
pl.colorbar() 

pl.subplot(223) 
pl.title('all data; max = {0:.1f}'.format(data.max())) 
pl.pcolormesh(lon, lat, np.transpose(data)) 
pl.colorbar() 

pl.subplot(224) 
pl.title('masked data; max = {0:.1f}'.format(data_m.max())) 
pl.pcolormesh(lon, lat, np.transpose(data_m)) 
pl.colorbar() 

was zur Folge hat: