-
Notifications
You must be signed in to change notification settings - Fork 0
/
Resampling_RACMO5km_to_CARRA_code.py
152 lines (111 loc) · 4.09 KB
/
Resampling_RACMO5km_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
147
148
149
150
151
152
# -*- 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
import gzip
import netCDF4
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 RACMO coordinates
#RACMO 5km data
fn=raw_path+'snowfall.1958-2020.FGRN055_BN_RACMO2.3p2_ERA5_3h_FGRN055.YY.nc.gz'
with gzip.open(fn) as gz:
ds=xr.open_dataset(gz,decode_times=False)
#load RACMO dataset
RACMO_data = np.array(ds.snowfall[0, 0, :, :])
lon=ds.lon.values
lat=ds.lat.values
# #from lat/lon to meters
# inProj = Proj(init='epsg:4326')
# outProj = Proj(init='epsg:3413')
niR=np.shape(lon)[0]
njR=np.shape(lon)[1]
x2, y2 = lon.flatten(), lat.flatten()
rx, ry = transform(inProj, outProj, x2, y2)
rx_mat = rx.reshape(niR, njR)
ry_mat = ry.reshape(niR, njR)
# x1, y1 = lon360_to_lon180(ds.x.values.flatten()), ds.y.values.flatten()
# rx, ry = transform(inProj, outProj, lon_mesh, lat_mesh)
# rx_mat = rx.reshape(niR, njR)
# ry_mat = ry.reshape(niR, njR)
cols2, rows2 = np.meshgrid(np.arange(njR), np.arange(niR))
# lat_m, lon_m = np.meshgrid(ds.LAT.values, ds.LON.values)
RACMO_positions = pd.DataFrame({
"row_r": rows2.ravel(),
"col_r": cols2.ravel(),
"lon_r": rx_mat.ravel(),
"lat_r": ry_mat.ravel(),})
# %% nearest neighbours
#dataset to be upscaled -> RACMO
nA = np.column_stack((rx_mat.ravel(), ry_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_RACMO = RACMO_positions.iloc[idx]
#%% Output resampling key ERA5 data in CARRA grid
CARRA_cells_for_RACMO.to_pickle(path_tools+'resampling_key_RACMO5km_to_CARRA.pkl')
#%% visualisation
new_grid= RACMO_data[CARRA_cells_for_RACMO.row_r, CARRA_cells_for_RACMO.col_r]
new_grid = new_grid.reshape(ni, nj)
plt.imshow(new_grid)
plt.colorbar()