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

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 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71
            if grid_type == 'z' or grid_type == 't':
                # extract example amplitude grid for Z, U and V and change NaNs to 0 (for land) and other values to 1 (for water)
                mask = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes'] + constituents[0] + '_Z.nc').variables['amplitude'][:])))
                mask[mask != 18446744073709551616.00000] = 1
                mask[mask == 18446744073709551616.00000] = 0
                self.mask_dataset[mz_name] = mask

                #read and convert the height_dataset file to complex and store in dicts
                hRe = []
                hIm = []
                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'][:])
                for ncon in range(len(constituents)):
                    amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+str(constituents[ncon])+'_Z.nc').variables['amplitude'][:])))
                    # set fill values to zero
                    amp[amp == 18446744073709551616.00000] = 0
72
                    # convert amp to m from cm
73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98
                    amp = amp/100.00
                    phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_Z.nc').variables['phase'][:])))
                    # set fill values to 0
                    phase[phase == 18446744073709551616.00000] = 0
                    # convert to real and imaginary conjugates and also convert to radians.
                    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]

            elif grid_type == 'u':
                mask = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes'] + constituents[0] + '_U.nc').variables['Ua'][:])))
                mask[mask != 18446744073709551616.00000] = 1
                mask[mask == 18446744073709551616.00000] = 0
                self.mask_dataset[mu_name] = mask
                #read and convert the velocity_dataset files to complex

                URe = []
                UIm = []
                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'][:])
                for ncon in range(len(constituents)):
                    amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ua'][:])))
                    # set fill values to zero
                    amp[amp == 18446744073709551616.00000] = 0
99
                    # convert amp units to m/s
100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
                    amp = amp/100.00
                    phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_U.nc').variables['Ug'][:])))
                    phase[phase == 18446744073709551616.00000] = 0
                    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]

            elif grid_type == 'v':
                mask = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes'] + constituents[0] + '_V.nc').variables['Va'][:])))
                mask[mask != 18446744073709551616.00000] = 1
                mask[mask == 18446744073709551616.00000] = 0
                self.mask_dataset[mv_name] = mask

                VRe = []
                VIm = []
                lat_v = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_V.nc').variables['lat'][:])
                lon_v = np.array(Dataset(settings['tide_fes'] + constituents[0] + '_V.nc').variables['lon'][:])
                for ncon in range(len(constituents)):
                    amp = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Va'][:])))
                    # set fill value to zero
                    amp[amp == 18446744073709551616.00000] = 0
123
                    # convert amp units to m/s
124 125 126 127 128 129 130 131
                    amp = amp/100.00
                    phase = np.ma.MaskedArray.filled(np.flipud(np.rot90(Dataset(settings['tide_fes']+constituents[ncon]+'_V.nc').variables['Vg'][:])))
                    phase[phase == 18446744073709551616.00000] = 0
                    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
132

thopri's avatar
thopri committed
133 134 135
            # 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']))
136 137
            # set fill values to zero
            height_z[height_z <0.00] = 0.00
thopri's avatar
thopri committed
138

thopri's avatar
thopri committed
139 140 141 142 143 144 145 146
        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
147
            lon_z = np.concatenate(([lon_z[0]-lon_resolution, ], lon_z,[lon_z[-1]+lon_resolution, ]))
thopri's avatar
thopri committed
148
            height_z = np.concatenate(([height_z[-1, :], ], height_z, [height_z[0, :],]), axis=0)
149
            mask_z = np.concatenate(([mask[-1, :], ], mask, [mask[0, :], ]), axis=0)
thopri's avatar
thopri committed
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 195 196 197

        #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
198
            self.amp, self.gph = self.interpolate_constituents(self.Uvelocity_dataset,
thopri's avatar
thopri committed
199 200 201
                                                               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
202
            self.amp, self.gph = self.interpolate_constituents(self.Vvelocity_dataset,
thopri's avatar
thopri committed
203 204 205 206 207 208 209 210 211
                                                               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
212 213
        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
214

215 216 217
        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
218
        #data = data.reshape(1,nc_dataset['M2'][real_var_name].shape)
thopri's avatar
thopri committed
219
        # data[data==0] = np.NaN
thopri's avatar
thopri committed
220

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

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

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

262 263 264 265
            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
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
        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)
324

thopri's avatar
thopri committed
325 326 327 328 329 330 331 332 333 334 335 336
    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')