fes_extract_HC.py 16.3 KB
Newer Older
thopri's avatar
thopri committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14
'''
This is to extract the tidal harmonic constants out of a tidal model
for a given locations
[amp,Gph] = tpxo_extract_HC(Model,lat,lon,type,Cid)

@author: Mr. Srikanth Nagella
'''

# pylint: disable=E1103
# pylint: disable=no-name-in-module
from netCDF4 import Dataset
from scipy import interpolate
import numpy as np

thopri's avatar
thopri committed
15
class HcExtract(object):
thopri's avatar
thopri committed
16
    """ This is FES model extract_hc.c implementation in python adapted from tpxo_extract_HC.py"""
thopri's avatar
thopri committed
17 18 19 20 21 22
    def __init__(self, settings, lat, lon, grid_type):
        """initialises the Extract of tide information from the netcdf
           Tidal files"""
	# Set tide model
        tide_model = 'fes'
        if tide_model == 'fes':  # Define stuff to generalise Tide model
thopri's avatar
thopri committed
23

24 25
            hRe_name = 'hRe'
            hIm_name = 'hIm'
thopri's avatar
thopri committed
26 27
            lon_z_name = 'lon_z'
            lat_z_name = 'lat_z'
28 29
            URe_name = 'URe'
            UIm_name = 'UIm'
thopri's avatar
thopri committed
30 31
            lon_u_name = 'lon_u'
            lat_u_name = 'lat_u'
32 33
            VRe_name = 'VRe'
            VIm_name = 'VIm'
thopri's avatar
thopri committed
34 35 36 37 38
            lon_v_name = 'lon_v'
            lat_v_name = 'lat_v'
            mz_name = 'mask_z'
            mu_name = 'mask_u'
            mv_name = 'mask_v'
thopri's avatar
thopri committed
39

40 41 42 43 44 45
            # create list of HC using namelist file as reference
            constituents = list(settings['clname'].values())
            # clean strings in list and change to upper case if not already
            for i in range(len(constituents)):
                constituents[i] = constituents[i].strip("',/\n")
                constituents[i] = constituents[i].upper()
thopri's avatar
thopri committed
46 47

            self.cons = constituents
thopri's avatar
thopri committed
48
            self.mask_dataset = {}
thopri's avatar
thopri committed
49 50

            # extract lon and lat z data
thopri's avatar
thopri committed
51 52
            lon_z = np.array(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['lon'])
            lat_z = np.array(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['lat'])
thopri's avatar
thopri committed
53 54 55
            lon_resolution = lon_z[1] - lon_z[0]
            data_in_km = 0 # added to maintain the reference to matlab tmd code

56
            # extract example amplitude grid for Z, U and V and change NaNs to 0 (for land) and other values to 1 (for water)
57
            mask_z = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_Z.nc').variables['amplitude'][:])))
58
            mask_z[mask_z != 18446744073709551616.00000] = 1
59
            mask_z[mask_z == 18446744073709551616.00000] = 0
thopri's avatar
thopri committed
60 61
            self.mask_dataset[mz_name] = mask_z

62
            mask_u = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_U.nc').variables['Ua'][:])))
63
            mask_z[mask_z != 18446744073709551616.00000] = 1
64
            mask_u[mask_u == 18446744073709551616.00000] = 0
thopri's avatar
thopri committed
65 66
            self.mask_dataset[mu_name] = mask_u

67
            mask_v = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[0]+'_V.nc').variables['Va'][:])))
68
            mask_z[mask_z != 18446744073709551616.00000] = 1
69
            mask_v[mask_v == 18446744073709551616.00000] = 0
thopri's avatar
thopri committed
70 71 72
            self.mask_dataset[mv_name] = mask_v

            #read and convert the height_dataset file to complex and store in dicts
73 74
            hRe = []
            hIm = []
thopri's avatar
thopri committed
75 76
            lat_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lat'][:])
            lon_z = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['lon'][:])
thopri's avatar
thopri committed
77
            for ncon in range(len(constituents)):
78
                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:])))
79
                # set fill values to zero
80
                amp[amp == 18446744073709551616.00000] = 0
81
                # convert to m from cm
82 83
                amp = amp/100.00
                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:])))
84
                phase[phase == 18446744073709551616.00000] = 0
85 86 87 88 89
                hRe.append(amp*np.cos(phase*(np.pi/180)))
                hIm.append(-amp*np.sin(phase*(np.pi/180)))
            hRe = np.stack(hRe)
            hIm = np.stack(hIm)
            self.height_dataset = [lon_z,lat_z,hRe,hIm]
thopri's avatar
thopri committed
90 91

            #read and convert the velocity_dataset files to complex
92 93
            URe = []
            UIm = []
thopri's avatar
thopri committed
94 95
            lat_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lat'][:])
            lon_u = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['lon'][:])
thopri's avatar
thopri committed
96
            for ncon in range(len(constituents)):
97
                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:])))
98
                # set fill values to zero
99
                amp[amp == 18446744073709551616.00000] = 0
100
                # convert to cm from m
101 102
                amp = amp/100.00
                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:])))
103
                phase[phase == 18446744073709551616.00000] = 0
104 105 106 107 108 109 110 111
                URe.append(amp*np.cos(phase*(np.pi/180)))
                UIm.append(-amp*np.sin(phase*(np.pi/180)))
            URe = np.stack(URe)
            UIm = np.stack(UIm)
            self.Uvelocity_dataset = [lon_u,lat_u,URe,UIm]

            VRe = []
            VIm = []
thopri's avatar
thopri committed
112 113
            lat_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_V.nc').variables['lat'][:])
            lon_v = np.array(Dataset(settings['tide_fes'] + constituents[ncon] + '_V.nc').variables['lon'][:])
thopri's avatar
thopri committed
114
            for ncon in range(len(constituents)):
115
                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:])))
116
                # set fill value to zero
117
                amp[amp == 18446744073709551616.00000] = 0
118
                # convert m to cm
119 120
                amp = amp/100.00
                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:])))
121
                phase[phase == 18446744073709551616.00000] = 0
122 123 124 125 126
                VRe.append(amp*np.cos(phase*(np.pi/180)))
                VIm.append(-amp*np.sin(phase*(np.pi/180)))
            VRe = np.stack(VRe)
            VIm = np.stack(VIm)
            self.Vvelocity_dataset = [lon_v,lat_v,VRe,VIm]
thopri's avatar
thopri committed
127

thopri's avatar
thopri committed
128 129 130
            # open grid variables these are resampled TPXO parameters so may not work correctly.
            self.grid = Dataset(settings['tide_fes']+'grid_fes.nc')
            height_z = np.array(np.rot90(self.grid.variables['hz']))
thopri's avatar
thopri committed
131

thopri's avatar
thopri committed
132 133 134 135 136 137 138 139
        else:
           print('Don''t know that tide model')

        # Wrap coordinates in longitude if the domain is global
        glob = 0
        if lon_z[-1]-lon_z[0] == 360-lon_resolution:
            glob = 1
        if glob == 1:
thopri's avatar
thopri committed
140
            lon_z = np.concatenate(([lon_z[0]-lon_resolution, ], lon_z,[lon_z[-1]+lon_resolution, ]))
thopri's avatar
thopri committed
141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190
            height_z = np.concatenate(([height_z[-1, :], ], height_z, [height_z[0, :],]), axis=0)
            mask_z = np.concatenate(([mask_z[-1, :], ], mask_z, [mask_z[0, :], ]), axis=0)

        #adjust lon convention
        xmin = np.min(lon)

        if data_in_km == 0:
            if xmin < lon_z[0]:
                lon[lon < 0] = lon[lon < 0] + 360
            if xmin > lon_z[-1]:
                lon[lon > 180] = lon[lon > 180]-360

        #height_z[height_z==0] = np.NaN
#        f=interpolate.RectBivariateSpline(lon_z,lat_z,height_z,kx=1,ky=1)
#        depth = np.zeros(lon.size)
#        for idx in range(lon.size):
#            depth[idx] = f(lon[idx],lat[idx])
#        print depth[369:371]

#        H2 = np.ravel(height_z)
#        H2[H2==0] = np.NaN
#        points= np.concatenate((np.ravel(self.height_dataset.variables['lon_z']),
#                                np.ravel(self.height_dataset.variables['lat_z'])))
#        points= np.reshape(points,(points.shape[0]/2,2),order='F')
#        print points.shape
#        print np.ravel(height_z).shape
#        depth = interpolate.griddata(points,H2,(lon,lat))
#        print depth
#        print depth.shape

        height_z[height_z == 0] = np.NaN
        lonlat = np.concatenate((lon, lat))
        lonlat = np.reshape(lonlat, (lon.size, 2), order='F')

        depth = interpolate.interpn((lon_z, lat_z), height_z, lonlat)
#        f=interpolate.RectBivariateSpline(lon_z,lat_z,mask_z,kx=1,ky=1)
#        depth_mask = np.zeros(lon.size)
#        for idx in range(lon.size):
#            depth_mask[idx] = f(lon[idx],lat[idx])
        depth_mask = interpolate.interpn((lon_z, lat_z), mask_z, lonlat)
        index = np.where((np.isnan(depth)) & (depth_mask > 0))

        if index[0].size != 0:
            depth[index] = bilinear_interpolation(lon_z, lat_z, height_z, lon[index], lat[index])

        if grid_type == 'z' or grid_type == 't':
            self.amp, self.gph = self.interpolate_constituents(self.height_dataset,
                                                               hRe_name, hIm_name, lon_z_name, lat_z_name,
							       lon, lat, maskname=mz_name)
        elif grid_type == 'u':
thopri's avatar
thopri committed
191
            self.amp, self.gph = self.interpolate_constituents(self.Uvelocity_dataset,
thopri's avatar
thopri committed
192 193 194
                                                               URe_name, UIm_name, lon_u_name, lat_u_name,
                                                               lon, lat, depth, maskname=mu_name)
        elif grid_type == 'v':
thopri's avatar
thopri committed
195
            self.amp, self.gph = self.interpolate_constituents(self.Vvelocity_dataset,
thopri's avatar
thopri committed
196 197 198 199 200 201 202 203 204
                                                               VRe_name, VIm_name, lon_v_name, lat_v_name,
                                                               lon, lat, depth, maskname=mv_name)
        else:
            print('Unknown grid_type')
            return

    def interpolate_constituents(self, nc_dataset, real_var_name, img_var_name, lon_var_name,
                                 lat_var_name, lon, lat, height_data=None, maskname=None):
        """ Interpolates the tidal constituents along the given lat lon coordinates """
thopri's avatar
thopri committed
205 206
        amp = np.zeros((len(nc_dataset[2]), lon.shape[0]))
        gph = np.zeros((len(nc_dataset[2]), lon.shape[0]))
thopri's avatar
thopri committed
207

208 209 210
        data = np.array(np.ravel(nc_dataset[2]), dtype=complex)
        data.imag = np.array(np.ravel(nc_dataset[3]))
        data = data.reshape(nc_dataset[2].shape)
thopri's avatar
thopri committed
211
        #data = data.reshape(1,nc_dataset['M2'][real_var_name].shape)
thopri's avatar
thopri committed
212
        # data[data==0] = np.NaN
thopri's avatar
thopri committed
213

thopri's avatar
thopri committed
214
        # Lat Lon values
thopri's avatar
thopri committed
215 216
        x_values = nc_dataset[0]
        y_values = nc_dataset[1]
thopri's avatar
thopri committed
217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235
        x_resolution = x_values[1] - x_values[0]
        glob = 0
        if x_values[-1]-x_values[0] == 360-x_resolution:
            glob = 1

        if glob == 1:
            x_values = np.concatenate(([x_values[0]-x_resolution,], x_values,
                                       [x_values[-1]+x_resolution, ]))

        #adjust lon convention
        xmin = np.min(lon)
        if xmin < x_values[0]:
            lon[lon < 0] = lon[lon < 0] + 360
        if xmin > x_values[-1]:
            lon[lon > 180] = lon[lon > 180]-360

        lonlat = np.concatenate((lon, lat))
        lonlat = np.reshape(lonlat, (lon.size, 2), order='F')

thopri's avatar
thopri committed
236
        mask = self.mask_dataset[maskname]
thopri's avatar
thopri committed
237 238 239 240
        mask = np.concatenate(([mask[-1, :], ], mask, [mask[0, :], ]), axis=0)
        #interpolate the mask values
        maskedpoints = interpolate.interpn((x_values, y_values), mask, lonlat)

241 242
        data_temp = np.zeros((data.shape[0], lon.shape[0], 2, ))
        for cons_index in range(data.shape[0]):
thopri's avatar
thopri committed
243
            #interpolate real values
244 245
            data_temp[cons_index, :, 0] = interpolate_data(x_values, y_values,
                                                                data[cons_index, :, :].real,
thopri's avatar
thopri committed
246 247
                                                                maskedpoints, lonlat)
            #interpolate imag values
248 249
            data_temp[cons_index, :, 1] = interpolate_data(x_values, y_values,
                                                                data[cons_index, :, :].imag,
thopri's avatar
thopri committed
250 251 252 253
                                                                maskedpoints, lonlat)

            #for velocity_dataset values
            if height_data is not None:
254 255
                data_temp[cons_index, :, 0] = data_temp[cons_index, :, 0]/height_data*100
                data_temp[cons_index, :, 1] = data_temp[cons_index, :, 1]/height_data*100
thopri's avatar
thopri committed
256

257 258
            zcomplex = np.array(data_temp[cons_index, :, 0], dtype=complex)
            zcomplex.imag = data_temp[cons_index, :, 1]
thopri's avatar
thopri committed
259

260 261 262 263
            amp[cons_index, :] = np.absolute(zcomplex)
            gph[cons_index, :] = np.arctan2(-1*zcomplex.imag, zcomplex.real)
        gph = gph*180.0/np.pi
        gph[gph < 0] = gph[gph < 0]+360.0
thopri's avatar
thopri committed
264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336
        return amp, gph

def interpolate_data(lon, lat, data, mask, lonlat):
    """ Interpolate data data on regular grid for given lonlat coordinates """
    result = np.zeros((lonlat.shape[0], ))
    data[data == 0] = np.NaN
    data = np.concatenate(([data[-1, :], ], data, [data[0, :], ]), axis=0)
    result[:] = interpolate.interpn((lon, lat), data, lonlat)
    index = np.where((np.isnan(result)) & (mask > 0))
    if index[0].size != 0:
        result[index] = bilinear_interpolation(lon, lat, data, np.ravel(lonlat[index, 0]),
                                               np.ravel(lonlat[index, 1]))
    return result

def bilinear_interpolation(lon, lat, data, lon_new, lat_new):
    """ Does a bilinear interpolation of grid where the data values are NaN's"""
    glob = 0
    lon_resolution = lon[1] - lon[0]
    if lon[-1] - lon[1] == 360 - lon_resolution:
        glob = 1
    inan = np.where(np.isnan(data))
    data[inan] = 0
    mask = np.zeros(data.shape)
    mask[data != 0] = 1
#    n = lon.size
#    m = lat.size
    if lon.size != data.shape[0] or lat.size != data.shape[1]:
        print('Check Dimensions')
        return np.NaN
    if glob == 1:
        lon = np.concatenate(([lon[0] - 2 * lon_resolution, lon[0] - lon_resolution, ],
                              lon, [lon[-1] + lon_resolution, lon[-1] + 2 * lon_resolution]))
        data = np.concatenate((data[-2, :], data[-1, :], data, data[0, :], data[1, :]), axis=0)
        mask = np.concatenate((mask[-2, :], mask[-1, :], mask, mask[0, :], mask[1, :]), axis=0)
    lon_new_copy = lon_new

    nonmask_index = np.where((lon_new_copy < lon[0]) & (lon_new_copy > lon[-1]))
    if lon[-1] > 180:
        lon_new_copy[nonmask_index] = lon_new_copy[nonmask_index] + 360
    if lon[-1] < 0:
        lon_new_copy[nonmask_index] = lon_new_copy[nonmask_index] - 360
    lon_new_copy[lon_new_copy > 360] = lon_new_copy[lon_new_copy > 360] - 360
    lon_new_copy[lon_new_copy < -180] = lon_new_copy[lon_new_copy < -180] + 360

    weight_factor_0 = 1 / (4 + 2 * np.sqrt(2))
    weight_factor_1 = weight_factor_0 / np.sqrt(2)
    height_temp = weight_factor_1 * data[0:-2, 0:-2] + weight_factor_0 * data[0:-2, 1:-1] + \
                  weight_factor_1 * data[0:-2, 2:] + weight_factor_1 * data[2:, 0:-2] + \
                  weight_factor_0 * data[2:, 1:-1] + weight_factor_1 * data[2:, 2:] + \
                  weight_factor_0 * data[1:-1, 0:-2] + weight_factor_0 * data[1:-1, 2:]
    mask_temp = weight_factor_1 * mask[0:-2, 0:-2] + weight_factor_0 * mask[0:-2, 1:-1] + \
                weight_factor_1 * mask[0:-2, 2:] + weight_factor_1 * mask[2:, 0:-2] + \
                weight_factor_0 * mask[2:, 1:-1] + weight_factor_1 * mask[2:, 2:] + \
                weight_factor_0 * mask[1:-1, 0:-2] + weight_factor_0 * mask[1:-1, 2:]
    mask_temp[mask_temp == 0] = 1
    data_copy = data.copy()
    data_copy[1:-1, 1:-1] = np.divide(height_temp, mask_temp)
    nonmask_index = np.where(mask == 1)
    data_copy[nonmask_index] = data[nonmask_index]
    data_copy[data_copy == 0] = np.NaN
    lonlat = np.concatenate((lon_new_copy, lat_new))
    lonlat = np.reshape(lonlat, (lon_new_copy.size, 2), order='F')
    result = interpolate.interpn((lon, lat), data_copy, lonlat)

    return result


#lat=[42.8920,42.9549,43.0178]
#lon=[339.4313,339.4324,339.4335]
#lat_u=[42.8916,42.9545,43.0174]
#lon_u=[339.4735,339.4746,339.4757]
#lat = np.array(lat_u)
#lon = np.array(lon_u)
337
#lon = TPXO_Extract(lat,lon,'velocity_dataset')