使用Python的NetCDF文件的循環緯度/經度作物

Question

我正在處理從國家氣象局ftp站點將原始二進制雷達數據導入到服務器的項目。使用天氣和氣候工具包數據導出工具，我將數據轉換爲netCDF文件。以下是對.NC文件「ncdump -h」命令的結果：

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]]" ; 
} 

我想找到的CREF變量，我可以與netCDF4和numpy的圖書館做很容易的最大入口在python：

import netCDF4 
import numpy 

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

不過，我希望能過濾的netCDF文件，以便最大和最小的僅在一個給定的經/緯度的一定距離發現。換句話說，我希望在指定的緯度/經度周圍「裁剪」一個指定半徑的圓。我發現how to crop a square通過另一個SO線程，但無法弄清楚一個圓圈是如何工作的。

Answer 1

我會計算中心和每個緯度/經度對（2D網格）之間的距離，並用它來構造一個可以應用於數據的蒙版。一旦屏蔽後，您可以再次使用numpy函數來計算統計信息，如max()。

例如，使用從https://stackoverflow.com/a/4913653/3581217的haversine()功能，修改爲矢量化版本，可以直接應用到numpy陣列：

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

導致：