swn_intrp.py

import numpy as np
from netCDF4 import Dataset
import sys

# create .grd file

xpc = 31.0
ypc = 71.0
xlenc = 18
ylenc = 9
mxc = 400
myc = 200
nx = mxc + 1
ny = myc + 1
shg_x = xlenc / mxc
shg_y = ylenc / myc

grd_x = np.zeros(nx)
grd_y = np.zeros(ny)
x = xpc
y = ypc

for i in range(int(nx)):
    grd_x[i] = x
    x += shg_x

for i in range(int(ny)):
    grd_y[i] = y
    y += shg_y

k = 0
grd = np.zeros((int(nx*ny), 2))
for i in range(int(nx)):
    for j in range(int(ny)):
        grd[k][0] = grd_x[i]
        grd[k][1] = grd_y[j]
        k += 1
np.savetxt('swn_grd2.dat', grd)

# grd_fname = 'swn_grd2.dat'
# grd_file = open(grd_fname, 'w')
# for i in range(nx):
#     for j in range(ny):
#         grd_file.write(str(grd_x[i])+'  '+str(grd_y[j])+'\n')

# open .nc file

longitude_name = 'lon'
latitude_name = 'lat'

nc_fileID2 = 'FED.nc'
root2 = Dataset(nc_fileID2)
u10_name = 'u10'
v10_name = 'v10'
u10v10_lat_name = 'latitude'
u10v10_lon_name = 'longitude'
long_u10v10 = root2.variables[u10v10_lon_name][:]
lat_u10v10 = root2.variables[u10v10_lat_name][:]

u10 = root2.variables[u10_name][:,:,:]
v10 = root2.variables[v10_name][:,:,:]

print(lat_u10v10)
print(long_u10v10)

w_o = np.zeros((nx*ny, 4))
u10v10_l = int(nx*ny)

time = 1

u10_swn = np.zeros((nx, ny))
v10_swn = np.zeros((nx, ny))

for t in range(time):
    fname_u10v10 = 'u10v10_'+str(t)+'.dat'
    u10v10_outf = open(fname_u10v10, 'w')
    for i in range(nx):
        for n in range(ny):
            ix = grd_x[i]
            iy = grd_y[n]
            for j in range(len(lat_u10v10)-1):
                for k in range(len(long_u10v10)-1):
                    if (long_u10v10[k] <= ix <= long_u10v10[k+1]):
                        if (lat_u10v10[j] >= iy >= lat_u10v10[j+1]):
                            Q_11_u = u10[t][j+1][k]
                            Q_12_u = u10[t][j][k]
                            Q_22_u = u10[t][j][k+1]
                            Q_21_u = u10[t][j+1][k+1]
                            x1, x2 = long_u10v10[k], long_u10v10[k+1]
                            y1, y2 = lat_u10v10[j+1], lat_u10v10[j]
                            block1 = (Q_11_u/((x2 - x1)*(y2 - y1)))*(x2 - ix)*(y2 - iy)
                            block2 = (Q_21_u/((x2 - x1)*(y2 - y1)))*(ix - x1)*(y2 - iy)
                            block3 = (Q_12_u/((x2 - x1)*(y2 - y1)))*(x2 - ix)*(iy - y1)
                            block4 = (Q_22_u/((x2 - x1)*(y2 - y1)))*(ix - x1)*(iy - y1)
                            f1 = block1 + block2 + block3 + block4
                            Q_11_v = v10[t][j+1][k]
                            Q_12_v = v10[t][j][k]
                            Q_22_v = v10[t][j][k+1]
                            Q_21_v = v10[t][j+1][k+1]
                            x1, x2 = long_u10v10[k], long_u10v10[k+1]
                            y1, y2 = lat_u10v10[j+1], lat_u10v10[j]
                            block1 = (Q_11_v/((x2 - x1)*(y2 - y1)))*(x2 - ix)*(y2 - iy)
                            block2 = (Q_21_v/((x2 - x1)*(y2 - y1)))*(ix - x1)*(y2 - iy)
                            block3 = (Q_12_v/((x2 - x1)*(y2 - y1)))*(x2 - ix)*(iy - y1)
                            block4 = (Q_22_v/((x2 - x1)*(y2 - y1)))*(ix - x1)*(iy - y1)
                            f2 = block1 + block2 + block3 + block4
            w_o[i][0] = ix
            w_o[i][1] = iy
            w_o[i][2] = f1
            w_o[i][3] = f2
            u10_swn[i][n] = f1
            v10_swn[i][n] = f2
            u10v10_outf.write(str(w_o[i][0])+' '+str(w_o[i][1])+' '+str(w_o[i][2])+' '+str(w_o[i][3])+'\n')
            print(i, n)
    w_swn_outname = 'w_swn_'+str(t)+'.dat'
    w_swnf = open(w_swn_outname, 'w')
    for q in range(nx):
        for w in range(ny):
            w_swnf.write(str(u10_swn[q][w])+' ')
        w_swnf.write('\n')
    for q in range(nx):
        for w in range(ny):
            w_swnf.write(str(v10_swn[q][w])+' ')
        w_swnf.write('\n')

# depth

longitude_name = 'lon'
latitude_name = 'lat'

nc_fileID = 'GRIDONE_2D_30.0_70.0_50.0_80.0.nc'
dpth_name = 'elevation'
root = Dataset(nc_fileID)

depth = root.variables[dpth_name][:, :]
long_depth = root.variables[longitude_name][:]
lat_depth = root.variables[latitude_name][:]

for i in range(depth.shape[0]):
    for j in range(depth.shape[1]):
        if (depth[i][j] < 0):
            depth[i][j] = -1 * depth[i][j]
        elif (depth[i][j] == 999):
            depth[i][j] = 1

depth_int = np.zeros((nx*ny, 3))
dpth_l = int(nx*ny) #  length of depth array

chk = input('interpolate depth? y/n')
if (chk == 'n'):
    sys.exit()

dp = np.zeros((nx, ny))
for i in range(nx):
    for n in range(ny):
        ix = grd_x[i]
        iy = grd_y[n]
        for j in range(len(lat_depth)-1):
            for k in range(len(long_depth)-1):
                if (long_depth[k] <= ix <= long_depth[k+1]):
                    if (lat_depth[j] <= iy <= lat_depth[j+1]):
                        Q_11 = depth[j+1][k]
                        Q_12 = depth[j][k]
                        Q_22 = depth[j][k+1]
                        Q_21 = depth[j+1][k+1]
                        x1, x2 = long_depth[k], long_depth[k+1]
                        y1, y2 = lat_depth[j+1], lat_depth[j]
                        block1 = (Q_11/((x2 - x1)*(y2 - y1)))*(x2 - ix)*(y2 - iy)
                        block2 = (Q_21/((x2 - x1)*(y2 - y1)))*(ix - x1)*(y2 - iy)
                        block3 = (Q_12/((x2 - x1)*(y2 - y1)))*(x2 - ix)*(iy - y1)
                        block4 = (Q_22/((x2 - x1)*(y2 - y1)))*(ix - x1)*(iy - y1)
                        f = block1 + block2 + block3 + block4
        depth_int[i][0] = ix
        depth_int[i][1] = iy
        depth_int[i][2] = f
        dp[i][n] = f
        print(i, n)
np.savetxt('dep_surf.dat', depth_int)
dep_swn_name = 'dep_swn.dat'
dep_swn = open(dep_swn_name, 'w')
for i in range(nx):
    for j in range(ny):
        dep_swn.write(str(dp[i][j])+' ')
    dep_swn.write('\n')