create_landseamask_MERRA.py

#!/usr/bin/env python

import os
import h5py
import numpy as np
import netCDF4 as nc

# ------------------------------ Download from NASA server and convert to hdf5 --------------------------------------------------
os.system('wget ftp://goldsmr2.sci.gsfc.nasa.gov/data/s4pa//MERRA_MONTHLY/MAC0NXASM.5.2.0/1979/MERRA300.prod.assim.const_2d_asm_Nx.00000000.hdf')
try:
        os.system('h4toh5 MERRA300.prod.assim.const_2d_asm_Nx.00000000.hdf')
except:
        exit('you need h4toh5 util, install h4h5tools from macports or equivalent')

# ------------------------------ Read fields we need from file --------------------------------------------------
f = h5py.File('MERRA300.prod.assim.const_2d_asm_Nx.00000000.h5','r')

frac_land    = f['EOSGRID']['Data Fields']['FRLAND'][:].squeeze()
frac_landice = f['EOSGRID']['Data Fields']['FRLANDICE'][:].squeeze()
frac_lake    = f['EOSGRID']['Data Fields']['FRLAKE'][:].squeeze()
frac_ocean   = f['EOSGRID']['Data Fields']['FROCEAN'][:].squeeze()

#area = f['EOSGRID']['Data Fields']['AREA'][:].squeeze()
lon  = f['EOSGRID']['Data Fields']['XDim'][:].squeeze()
lat  = f['EOSGRID']['Data Fields']['YDim'][:].squeeze()

f.close()

# ------------------------------ Create our binary land sea mask --------------------------------------------------
ny, nx = frac_ocean.shape
lsm = np.zeros((ny,nx))

threshold = 0.95
lsm[np.where(frac_ocean >= threshold)] = 1

# ------------------------------ Move to 0-360 --------------------------------------------------------------------
# find llon = 0
indx = np.abs(lon).argmin()

def offset_1d(a_in,index,lon=False):
        sx = a_in.shape[0]
        tx = sx - index
        a_out = a_in.copy()
        a_out[:] = -9999.
        a_out[:tx] = a_in[index:]
        if lon is True:
                a_out[tx:] = a_in[0:index] + 360.
        else:
                a_out[tx:] = a_in[0:index]
        return a_out

def offset_2d(a_in,index):
        sx = a_in.shape[1]
        tx = sx - index
        a_out = a_in.copy()
        a_out[:,:] = -9999.
        a_out[:,:tx] = a_in[:,index:]
        a_out[:,tx:] = a_in[:,0:index]
        return a_out

lonout = offset_1d(lon,indx,lon=True)
latout = lat.copy()

lsmout = offset_2d(lsm,indx)

# ------------------------------ Write the mask to netcdf ---------------------------------------------------------
fout = nc.Dataset('./MERRA_lsm.nc','w',format='NETCDF3_CLASSIC')
fout.description = 'MERRA land sea mask '
# dimensions
fout.createDimension('lat', ny)
fout.createDimension('lon', nx)
# variables
latitudes  = fout.createVariable('lat', 'f8', ('lat',))
longitudes = fout.createVariable('lon', 'f8', ('lon',))
mask       = fout.createVariable('lsm', 'i4', ('lat','lon',))
# data
latitudes[:]    = latout
longitudes[:]   = lonout
mask[:,:]       = lsmout
# close
fout.close()