fes_extract_HC.py 16.5 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

thopri's avatar
thopri committed
24 25
            hRe_name = 'amp_stack'
            hIm_name = 'pha_stack'
thopri's avatar
thopri committed
26 27
            lon_z_name = 'lon_z'
            lat_z_name = 'lat_z'
thopri's avatar
thopri committed
28 29
            URe_name = 'amp_stack'
            UIm_name = 'pha_stack'
thopri's avatar
thopri committed
30 31
            lon_u_name = 'lon_u'
            lat_u_name = 'lat_u'
thopri's avatar
thopri committed
32 33
            VRe_name = 'amp_stack'
            VIm_name = 'pha_stack'
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'][:])))
thopri's avatar
thopri committed
63
            mask_u[mask_u != 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'][:])))
thopri's avatar
thopri committed
68
            mask_v[mask_v != 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
thopri's avatar
thopri committed
73 74
            amp_stack = []
            pha_stack = []
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
thopri's avatar
thopri committed
85 86 87 88 89 90 91
                # convert from -180-180 to 0-360
                phase[phase < 0.0] = phase[phase < 0.0] + 360.0
                amp_stack.append(amp)
                pha_stack.append(phase)
            amp_stack = np.stack(amp_stack)
            pha_stack = np.stack(pha_stack)
            self.height_dataset = [lon_z,lat_z,amp_stack,pha_stack]
thopri's avatar
thopri committed
92 93

            #read and convert the velocity_dataset files to complex
thopri's avatar
thopri committed
94 95
            amp_stack = []
            pha_stack = []
thopri's avatar
thopri committed
96 97
            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
98
            for ncon in range(len(constituents)):
99
                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:])))
100
                # set fill values to zero
101
                amp[amp == 18446744073709551616.00000] = 0
102
                # convert to cm from m
103 104
                amp = amp/100.00
                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:])))
105
                phase[phase == 18446744073709551616.00000] = 0
thopri's avatar
thopri committed
106 107 108 109 110 111 112 113 114
                phase[phase < 0.0] = phase[phase < 0.0] + 360.0
                amp_stack.append(amp)
                pha_stack.append(phase)
            amp_stack = np.stack(amp_stack)
            pha_stack = np.stack(pha_stack)
            self.Uvelocity_dataset = [lon_u,lat_u,amp_stack,pha_stack]

            amp_stack = []
            pha_stack = []
thopri's avatar
thopri committed
115 116
            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
117
            for ncon in range(len(constituents)):
118
                amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:])))
119
                # set fill value to zero
120
                amp[amp == 18446744073709551616.00000] = 0
121
                # convert m to cm
122 123
                amp = amp/100.00
                phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:])))
124
                phase[phase == 18446744073709551616.00000] = 0
thopri's avatar
thopri committed
125 126 127 128 129 130
                phase[phase < 0.0] = phase[phase < 0.0] + 360.0
                amp_stack.append(amp)
                pha_stack.append(phase)
            amp_stack = np.stack(amp_stack)
            pha_stack = np.stack(pha_stack)
            self.Vvelocity_dataset = [lon_v,lat_v,amp_stack,pha_stack]
thopri's avatar
thopri committed
131

thopri's avatar
thopri committed
132 133 134
            # 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
135

thopri's avatar
thopri committed
136 137 138 139 140 141 142 143
        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
144
            lon_z = np.concatenate(([lon_z[0]-lon_resolution, ], lon_z,[lon_z[-1]+lon_resolution, ]))
thopri's avatar
thopri committed
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 191 192 193 194
            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
195
            self.amp, self.gph = self.interpolate_constituents(self.Uvelocity_dataset,
thopri's avatar
thopri committed
196 197 198
                                                               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
199
            self.amp, self.gph = self.interpolate_constituents(self.Vvelocity_dataset,
thopri's avatar
thopri committed
200 201 202 203 204 205 206 207 208
                                                               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
209 210
        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
211

thopri's avatar
thopri committed
212 213 214 215
        data_amp = np.array(np.ravel(nc_dataset[2]))
        data_pha = np.array(np.ravel(nc_dataset[3]))
        data_amp = data_amp.reshape(nc_dataset[2].shape)
        data_pha = data_pha.reshape(nc_dataset[2].shape)
thopri's avatar
thopri committed
216
        #data = data.reshape(1,nc_dataset['M2'][real_var_name].shape)
thopri's avatar
thopri committed
217
        # data[data==0] = np.NaN
thopri's avatar
thopri committed
218

thopri's avatar
thopri committed
219
        # Lat Lon values
thopri's avatar
thopri committed
220 221
        x_values = nc_dataset[0]
        y_values = nc_dataset[1]
thopri's avatar
thopri committed
222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240
        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
241
        mask = self.mask_dataset[maskname]
thopri's avatar
thopri committed
242 243 244 245
        mask = np.concatenate(([mask[-1, :], ], mask, [mask[0, :], ]), axis=0)
        #interpolate the mask values
        maskedpoints = interpolate.interpn((x_values, y_values), mask, lonlat)

thopri's avatar
thopri committed
246 247 248
        data_temp_amp = np.zeros((data_amp.shape[0], lon.shape[0], 1, ))
        data_temp_pha = np.zeros((data_pha.shape[0], lon.shape[0], 1,))
        for cons_index in range(data_amp.shape[0]):
thopri's avatar
thopri committed
249
            #interpolate real values
thopri's avatar
thopri committed
250 251
            data_temp_amp[cons_index, :, 0] = interpolate_data(x_values, y_values,
                                                                data_amp[cons_index, :, :],
thopri's avatar
thopri committed
252 253
                                                                maskedpoints, lonlat)
            #interpolate imag values
thopri's avatar
thopri committed
254 255
            data_temp_pha[cons_index, :, 0] = interpolate_data(x_values, y_values,
                                                                data_pha[cons_index, :, :],
thopri's avatar
thopri committed
256 257 258 259
                                                                maskedpoints, lonlat)

            #for velocity_dataset values
            if height_data is not None:
thopri's avatar
thopri committed
260 261
                data_temp_amp[cons_index, :, 0] = data_temp_amp[cons_index, :, 0]/height_data*100
                data_temp_pha[cons_index, :, 0] = data_temp_pha[cons_index, :, 0]/height_data*100
thopri's avatar
thopri committed
262

thopri's avatar
thopri committed
263 264
            amp[cons_index, :] = data_temp_amp[cons_index, :, 0]
            gph[cons_index, :] = data_temp_pha[cons_index, :, 0]
thopri's avatar
thopri committed
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 337 338 339

        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)
#lon = TPXO_Extract(lat,lon,'velocity_dataset')