-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathResampling_JRA_to_CARRA_code.py
146 lines (105 loc) · 4 KB
/
Resampling_JRA_to_CARRA_code.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
# -*- coding: utf-8 -*-
"""
@author: Adrien Wehrlé, GEUS (Geological Survey of Denmark and Greenland)
Few modifications needed to be run with other datasets
-> modified by Armin Dachauer
"""
import numpy as np
import pandas as pd
import os
from pyproj import Proj, transform
import xarray as xr
from scipy.spatial import cKDTree
import matplotlib.pyplot as plt
AW = 0
AD=1
base_path = '/Users/jason/Dropbox/CARRA/CARRA_rain/'
if AW:
base_path = 'C:/Users/Pascal/Desktop/GEUS_2019/SICE_AW_JEB/SICE_AW_JEB/'\
+ 'CARRA_rain/'
if AD:
base_path='C:/Users/Armin/Documents/Work/GEUS/Github/CARRA_rain/'
raw_path='C:/Users/Armin/Documents/Work/GEUS/Github/CARRA/'
path_tools='C:/Users/Armin/Documents/Work/GEUS/Github/CARRA_tools/'
os.chdir(base_path)
# %% CARRA coordinates
def lon360_to_lon180(lon360):
# reduce the angle
lon180 = lon360 % 360
# force it to be the positive remainder, so that 0 <= angle < 360
lon180 = (lon180 + 360) % 360;
# force into the minimum absolute value residue class, so that -180 < angle <= 180
lon180[lon180 > 180] -= 360
return lon180
# CARRA West grid dims
ni = 1269 ; nj = 1069
# read lat lon arrays
fn = './ancil/2.5km_CARRA_west_lat_1269x1069.npy'
lat = np.fromfile(fn, dtype=np.float32)
clat_mat = lat.reshape(ni, nj)
fn = './ancil/2.5km_CARRA_west_lon_1269x1069.npy'
lon = np.fromfile(fn, dtype=np.float32)
lon_pn = lon360_to_lon180(lon)
clon_mat = lon_pn.reshape(ni, nj)
fn='./ancil/CARRA_W_domain_ice_mask.nc'
ds=xr.open_dataset(fn)
mask = np.array(ds.z)
# %% reproject 4326 (lat/lon) CARRA coordinates to 3413 (orth polar projection in meters)
#from lat/lon to meters
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:3413')
x1, y1 = clon_mat.flatten(), clat_mat.flatten()
cx, cy = transform(inProj, outProj, x1, y1)
cx_mat = cx.reshape(ni, nj)
cy_mat = cy.reshape(ni, nj)
cols, rows = np.meshgrid(np.arange(np.shape(clat_mat)[1]),
np.arange(np.shape(clat_mat)[0]))
CARRA_positions = pd.DataFrame({'rowc': rows.ravel(),
'colc': cols.ravel(),
'xc': cx_mat.ravel(),
'yc': cy_mat.ravel(),
'maskc': mask.flatten()})
#import CARRA datset
ds = xr.open_dataset(raw_path+'tp_2012.nc')
CARRA_data = np.array(ds.tp[0, :, :]).flatten()
# %% load JRA coordinates
#JRA data
fn=raw_path+'jra55snow9820.nc'
ds=xr.open_dataset(fn)
#from lat/lon to meters
inProj = Proj(init='epsg:4326')
outProj = Proj(init='epsg:3413')
niJ=np.shape(ds.lon.values)[0]
njJ=np.shape(ds.lat.values)[0]
lat_mesh, lon_mesh = np.meshgrid(ds.lat.values, ds.lon.values)
# x1, y1 = lon360_to_lon180(ds.lon.values.flatten()), ds.lat.values.flatten()
x1, y1 = lon360_to_lon180(lon_mesh.flatten()), lat_mesh.flatten()
jx,jy = transform(inProj, outProj, x1, y1)
jx_mat = jx.reshape(niJ, njJ)
jy_mat = jy.reshape(niJ, njJ)
# cols2, rows2 = np.meshgrid(np.arange(niJ), np.arange(njJ))
# lat_m, lon_m = np.meshgrid(ds.LAT.values, ds.LON.values)
cols2, rows2 = np.meshgrid(np.arange(np.shape(ds.lat.values)[0]), np.arange(np.shape(ds.lon.values)[0]))
lat_j, lon_j = np.meshgrid(ds.lat.values, ds.lon.values)
JRA_positions = pd.DataFrame({
"row_j": rows2.ravel(),
"col_j": cols2.ravel(),
"lon_j": lat_j.ravel(),
"lat_j": lon_j.ravel(),})
#load JRA dataset
JRA_data = np.array(ds.snow[:, :])
# %% nearest neighbours
#dataset to be upscaled -> JRA
nA = np.column_stack((jx_mat.ravel(), jy_mat.ravel()) )
#dataset to provide the desired grid -> CARRA
nB = np.column_stack((cx_mat.ravel(), cy_mat.ravel()))
btree = cKDTree(nA) #train dataset
dist, idx = btree.query(nB, k=1) #apply on grid
#collocate matching cells
CARRA_cells_for_JRA = JRA_positions.iloc[idx]
#%% Output resampling key ERA5 data in CARRA grid
CARRA_cells_for_JRA.to_pickle(path_tools+'resampling_key_JRA_to_CARRA.pkl')
#%% visualisation
new_grid= JRA_data[CARRA_cells_for_JRA.col_j, CARRA_cells_for_JRA.row_j]
new_grid = new_grid.reshape(ni, nj)
plt.imshow(new_grid)