Круговая лат / лон обрезка файла NetCDF с Python

Question

Я работаю над проектом, который импортирует необработанные двоичные радарные данные с ftp-сайта Национальной метеорологической службы на сервер. Используя инструмент экспорта данных Weather and Climate Toolkit, я конвертирую данные в файл netCDF. Ниже приведен результат выполнения команды "ncdump -h" в файле.nc:

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, чтобы максимальные и минимальные значения находились только на определенном расстоянии от заданного широты / долготы. Другими словами, я надеюсь "обрезать" круг указанного радиуса вокруг указанного широты / долготы. Я нашел, как обрезать квадрат через другой поток SO, но не могу понять, как будет работать круг.

Answer 1

Я бы вычислил расстояние между центром и каждой парой широта / долгота (двумерная сетка) и использовал это для построения маски, которую вы можете применить к своим данным. После маски можно снова использовать numpy функции для расчета статистики, как max(),

Например, используя haversine() функция из /questions/46704283/formula-haversine-v-python-azimut-i-rasstoyanie-mezhdu-dvumya-tochkami-gps/46704303#46704303, измененная в векторную версию, которую вы можете напрямую применить к 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()

Что приводит к: