interpolation1.py

import sys
import numpy as np
from netCDF4 import Dataset

# ......................................................................
# STEP 1
# read .grd file and take coordinates of grid points
# ......................................................................

# print grid
# input grid filename
print('\nstart read grid coordinates\n')
# grd_fname = input('.grid file name (without format): ') + '.grd'
grd_fname = 'Grid_006_.grd'
try:
    file_chk = open(grd_fname)
except FileNotFoundError:
    sys.exit('\nERROR! file not found')

# read .grd file and parameters
grd_s = []
file = open(grd_fname)
for line in file:
    grd_s.append(line)
crd_inf = grd_s[6].split()
lng = int(crd_inf[0])
lt = int(crd_inf[1])
head_L = 8
for i in range(head_L):
    del grd_s[0]
m = lng//5
if (lng%5 != 0):
    m += 1

# write array
x = []
x_app = []
y = []
y_app = []
stps_get = m*lt

for i in range(stps_get):
    if (i%m == 0):
        cx = grd_s[i].split()
        cy = grd_s[i + stps_get].split()
        for j in range(len(cx)):
            x_app.append(cx[j])
        for j in range(len(cy)):
            y_app.append(cy[j])
        for k in range(2):
            del x_app[0]
            del y_app[0]
        for j in range(len(x_app)):
            x.append(x_app[j])
        for j in range(len(y_app)):
            y.append(y_app[j])
        x_app.clear()
        y_app.clear()
    else:
        cx = grd_s[i].split()
        cy = grd_s[i + stps_get].split()
        for j in range(len(cx)):
            x_app.append(cx[j])
        for j in range(len(cy)):
            y_app.append(cy[j])
        for j in range(len(x_app)):
            x.append(x_app[j])
        for j in range(len(y_app)):
            y.append(y_app[j])
        x_app.clear()
        y_app.clear()
x_crd_d3dgrd = np.zeros(len(x))
y_crd_d3dgrd = np.zeros(len(y))
for i in range(len(x)):
    x_crd_d3dgrd[i] = float(x[i])
for i in range(len(y)):
    y_crd_d3dgrd[i] = float(y[i])
fgrid_out_name = 'grd_out1.dat'
fgrid_out = open(fgrid_out_name, 'w')

# write in file
# for i in range(len(x_crd_d3dgrd)):
#     fgrid_out.write(str(x_crd_d3dgrd[i]) + ' ' + str(y_crd_d3dgrd[i]) + '\n')
print('\nread .grd complete')

# ......................................................................
# STEP 2
# open and read .nc file
# read variables
# DOWNLOAD NETCDF4 DATA FILE FROM:
# http://apps.ecmwf.int/datasets/data/interim-full-daily/levtype=sfc/
# This script reads netCDF files. NetCDF (Network Common Data Form) is a set of software
# libraries and self-describing, machine-independent data formats that support the creation,
# access, and sharing of array-oriented scientific data.
# ......................................................................

print('\nstart read .nc file\n')

# sys.stdout.write('-'+key+' ')
# nc_fileID = str(input('.nc file name (without format): ')) + '.nc'
nc_fileID = 'dat_3.nc'
try:
    file_chk = open(nc_fileID)
except FileNotFoundError:
    sys.exit('\nERROR! file not found')

root = Dataset(nc_fileID)
dims = root.dimensions
ndims = len(dims)

dic = {}
dic2 = {}
print('\navailable variables in selected .NC file:\n')
vars = root.variables
nvars = len(vars)
n = 0
for var in vars:
    # sys.stdout.write('-'+var+' ')
    print('#',n,'   ',var, vars[var].shape)
    dic[str(var)] = n
    l = vars[var].shape
    dic2[str(var)] = len(l)
    n += 1
print('\n')

nc_data = []
iter_var = []
iter_var.append(input('write the variables you want to read (through the gap): ').split(' '))
iter_var = iter_var[0]
try:
    for v in iter_var:
        dic[v]
except KeyError:
    sys.exit('\nvariable name error\n')

var_inter_n = {}
ni = 0
for i in iter_var:
    if (dic2[i] == 1):
        nc_data.append(np.array(root.variables[i][:], dtype=np.float32))
    if (dic2[i] == 2):
        nc_data.append(np.array(root.variables[i][:,:], dtype=np.float32))
    if (dic2[i] == 3):
        nc_data.append(np.array(root.variables[i][:,:,:], dtype=np.float32))
    var_inter_n[i] = ni
    ni += 1

print('\nread complete\n')

# ......................................................................
# STEP 3
# Bilinear interpolation
# ......................................................................

# longitude_name = str(input('longitude variable name: '))
# latitude_name = str(input('latitude variable name: '))
# time_name = str(input('time variable name: '))
longitude_name = 'longitude'
latitude_name = 'latitude'
time_name = 'time'

# check longitude and latitude
m = 0  # latitude
n = 0  # longitude
ex_lt = 1
ex_lng = 1
for i in iter_var:
    if (i == latitude_name):
        ex_lt = 0
        break
    m += 1
for i in iter_var:
    if (i == longitude_name):
        ex_lng = 0
        break
    n += 1
if (ex_lng == 1 and ex_lt == 1):
    sys.exit('\nno longitude and latitude data, interpolation is impossible')
if (ex_lng == 1):
    sys.exit('\nno longitude data, interpolation is impossible')
if (ex_lt == 1):
    sys.exit('\nno latitude data, interpolation is impossible')
iter_var2 = []
for i in iter_var:
    if (i != longitude_name and i != latitude_name and i != time_name):
        iter_var2.append(i)
if (len(iter_var2) == 0):
    sys.exit('\nno data for interpolation')

# WARNING
# iter_var2 - readed variables without long, lat and time
# var_inter_n - number of interpolated variable (without lon, lat, time) in nc_data
# x\y_crd_d3dgrd - coordiantes of delft grid
# lat\long - coordinates of netcdf file grid
# .......

long = nc_data[n]
lat = nc_data[m]

for i in range(5):
    ix = x_crd_d3dgrd[i]
    iy = y_crd_d3dgrd[i]
    for j in range(len(long)-1):
        if (long[j] < ix < long[j+1]):
            # print(long[j], long[j+1], ix, j, j-1)
            for k in range(len(lat)-1):
                if (lat[k+1] < iy < lat[k]):
                    # print(long[j], long[j+1], lat[k], lat[k+1], ix, iy, j, k)
                    continue

int_u10 = np.array(len(x_crd_d3dgrd))
int_v10 = np.array(len(y_crd_d3dgrd))

# print(nc_data[0])
# print(nc_data[1])

if (len(x_crd_d3dgrd) == len(y_crd_d3dgrd)):
    crd_d3d = x_crd_d3dgrd
    l_crd_d3d = int(len(crd_d3d))
else:
    sys.exit('\ncrd d3d error')

f_int_all = np.zeros((l_crd_d3d, 3))

time = input('enter time: ')

for v in iter_var2:
    # for t in range(len(time)):
    for t in range(int(time)):
        for i in range(l_crd_d3d):
        # for i in len(crd_d3d):
            ix = x_crd_d3dgrd[i]
            iy = y_crd_d3dgrd[i]
            for j in range(len(long)-1):
                if (ix < -999 or iy < -999):
                    f_int = -999.999
                else:
                    if (long[j] <= ix <= long[j+1]):
                        # print(long[j], long[j+1], ix, j, j-1)
                        for k in range(len(lat)-1):
                            if (lat[k] >= iy >= lat[k+1]):
                                # nc_data[# variable][time, longitude, latitude]
                                # k - lat - y +
                                # j - long - x
                                # print(long[j], long[j+1], lat[k], lat[k+1], ix, iy, j, k)
                                # print(ix, iy, j, k)
                                Q_11 = nc_data[var_inter_n[v]][t, k, j]
                                Q_12 = nc_data[var_inter_n[v]][t, k+1, j]
                                Q_21 = nc_data[var_inter_n[v]][t, k, j+1]
                                Q_22 = nc_data[var_inter_n[v]][t, k+1, j+1]
                                x1, x2, y1, y2 = long[j], long[j+1], lat[k], lat[k+1]
                                block1 = (Q_11 * (x2 - ix) * (y2 - iy)) / ((x2 - x1) * (y2 - y1))
                                block2 = (Q_21 * (ix - x1) * (y2 - iy)) / ((x2 - x1) * (y2 - y1))
                                block3 = (Q_12 * (x2 - ix) * (iy - y1)) / ((x2 - x1) * (y2 - y1))
                                block4 = (Q_22 * (ix - x1) * (iy - y1)) / ((x2 - x1) * (y2 - y1))
                                f_int = block1 + block2 + block3 + block4
            f_int_all[i][0] = ix
            f_int_all[i][1] = iy
            f_int_all[i][2] = f_int
        fname_f = str(v) + '_' + str(t) + '.dat'
        # for i in range(l_crd_d3d):

        file1 = open(fname_f, 'w')
        for i in range(int(l_crd_d3d)):
            file1.write(str(f_int_all[i][0]) + ' ' + str(f_int_all[i][1]) + ' ' + str(f_int_all[i][2]) + '\n')

        # np.savetxt(fname_f, f_int_all, delimiter='  ', fmt='%1.5f %1.5f %1.5f')
        # print(f_int_all)
    # print(v, nc_data[var_inter_n[v]][1,1,1])

# longitude latitude time u10 v10 sp

# print(nc_data[3][1,:,:])
# print('_______')
# print(nc_data[3][1,0,0])
# print(nc_data[3][1,0,1])
# print(nc_data[3][1,0,2])