Demo Python script to write and read netCDF4 files

readWriteNetCDF.py

"""
DEMO TO WRITE AND READ NETCDF FILES
version 27 March 2018
by Antti Lipponen

Copyright (c) 2018 Antti Lipponen

Permission is hereby granted, free of charge, to any person obtaining
a copy of this software and associated documentation files (the
"Software"), to deal in the Software without restriction, including
without limitation the rights to use, copy, modify, merge, publish,
distribute, sublicense, and/or sell copies of the Software, and to
permit persons to whom the Software is furnished to do so, subject to
the following conditions:
The above copyright notice and this permission notice shall be included
in all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
"""

from netCDF4 import Dataset, num2date, date2num
from datetime import datetime
import numpy as np
import os

# Aim is to write & read netCDF4 file which contains distances from different towns
# 1 Jan 2018: the distance is from Kuopio, Finland
# 2 Jan 2018: the distance is from Sydney, Australia

# our regular lat/lon grid
lat = np.linspace(-90, 90, 180)
lon = np.linspace(-180, 180, 360)

# times
time = [datetime(2018, 1, 1), datetime(2018, 1, 2)]


# Function to compute distances
def dist(lat1, lon1, lat2, lon2):
    lat1R, lat2R, lon1R, lon2R = np.deg2rad(lat1), np.deg2rad(lat2), np.deg2rad(lon1), np.deg2rad(lon2)
    dlon = lon2R - lon1R
    dlat = lat2R - lat1R
    R = 6378.1
    a = (np.sin(dlat / 2.0))**2 + np.cos(lat1R) * np.cos(lat2R) * (np.sin(dlon / 2.0))**2
    c = 2.0 * np.arctan2(np.sqrt(a), np.sqrt(1.0 - a))
    d = R * c
    return d


# make grid
gridlon, gridlat = np.meshgrid(lon, lat)

# compute distances
lat_Kuopio, lon_Kuopio = 62.89, 27.68
distanceFromKuopio = dist(lat_Kuopio, lon_Kuopio, gridlat.ravel(), gridlon.ravel()).reshape(gridlon.shape)
lat_Sydney, lon_Sydney = -33.86, 151.22
distanceFromSydney = dist(lat_Sydney, lon_Sydney, gridlat.ravel(), gridlon.ravel()).reshape(gridlon.shape)

distanceData = np.stack((distanceFromKuopio, distanceFromSydney), axis=0)  # dimension (time, lat, lon) = (2, 180, 360)

# write the netCDF file
ncout = Dataset('distances.nc', 'w', format='NETCDF4')

# Add some attributes
ncout.History = 'File generated on {} (UTC) by {}'.format(datetime.utcnow().strftime('%c'), os.path.basename(__file__))
ncout.Description = 'Demo file containing distance data from Kuopio, Finland & Sydney, Australia'

# create dimensions
ncout.createDimension('time', len(time))
ncout.createDimension('latitude', len(lat))
ncout.createDimension('longitude', len(lon))

# save coordinates
ncout_latitude = ncout.createVariable('latitude', 'f4', ('latitude',))
ncout_latitude[:] = lat
ncout_longitude = ncout.createVariable('longitude', 'f4', ('longitude',))
ncout_longitude[:] = lon

# create & save time
units = 'seconds since 1980-01-01'
ncout_time = ncout.createVariable('time', 'i4', ('time',))
ncout_time[:] = date2num(time, units)
ncout_time.units = units

# save data
ncout_dist = ncout.createVariable('distance', 'f4', ('time', 'latitude', 'longitude'))
ncout_dist.long_name = 'Distance'
ncout_dist.units = 'km'
ncout_dist[:] = distanceData

ncout.close()

# Reading netcdf
ncin = Dataset('distances.nc', 'r')

# get description
description = ncin.Description  # or getattr(ncin, 'Description')
print('Description: {}'.format(description))

ncin_lat, ncin_lon = ncin['latitude'][:], ncin['longitude'][:]
ncin_distancedata = ncin['distance'][:]

ncin_time = num2date(ncin['time'][:], ncin['time'].units)

ncin.close()

print('Time instants in NC-file:')
for t in ncin_time:
        print('  {}'.format(t))

print('Data array dimension: {}'.format(ncin_distancedata.shape))