Circular lat/lon crop of a NetCDF file with Python

Question

I'm working on a project which imports raw binary radar data from the National Weather Service ftp site to a server. Using an the Weather and Climate Toolkit data-export tool, I convert the data into a netCDF file. The following is the result of a "ncdump -h" command on the .nc file: ``` 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…

Accepted Answer

I would calculate the distance between the center and each lat/lon pair (2D grid), and use that to construct a mask which you can apply to your data. Once masked, you can again simply use `numpy` functions to calculate statistics like `max()`.

For example, using the `haversine()` function from https://stackoverflow.com/a/4913653/3581217, modified to a vectorized version which you can directly apply onto `numpy` arrays:

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

Which results in:

Circular lat/lon crop of a NetCDF file with Python

Solution

Problem

Related problems