rho_indian_ocean.py

import xarray
import numpy
from dask.diagnostics import ProgressBar

in_path = '/g/data/w35/saw562/b.e11.B1850C5CN.f19_g16.0850cntl.001.pop.h.RHO.185001-200512.nc'
out_path = '/g/data/w35/saw562/test.nc'

min_lat = -35.5
max_lat = 0.5
min_lon = 18.5
max_lon = 130.5

# Chunks are found from the file using `ncdump -hs`:
# 		RHO:_ChunkSizes = 1, 30, 192, 160 ;
ds = xarray.open_dataset(in_path, chunks={'time': 1, 'z_t': 30, 'nlat': 192, 'nlon': 160})

# Lat and Lon are 2d - create a mask of the area we're interested in, then get the array indexes from that mask
lat_mask = numpy.logical_and(min_lat < ds.TLAT, ds.TLAT < max_lat)
lon_mask = numpy.logical_and(min_lon < ds.TLONG, ds.TLONG < max_lon)
lat_indices, lon_indices = numpy.where(numpy.logical_and(lat_mask, lon_mask))

# Use the min and max of the mask indices to select the region of interest
ds = ds.isel(nlat = slice(lat_indices.min(), lat_indices.max()),
             nlon = slice(lon_indices.min(), lon_indices.max()))

# Save the variables of interest, using compression
with ProgressBar():
  ds[['RHO','dz']].to_netcdf(out_path, 
                   encoding={'RHO': {
                     'zlib': True, 'shuffle': True, 'complevel': 4, 
                     'chunksizes': (1, 30, ds.RHO.shape[2], ds.RHO.shape[3])
                   }})