Gage adjustment¶

In [1]:

import os
import math
import copy
import warnings

import numpy as np
import numpy.ma as ma

import matplotlib.pyplot as plt
import cartopy.crs as ccrs
import xarray as xr

import wradlib as wrl

try:
    get_ipython().run_line_magic("matplotlib inline")
except:
    plt.ion()

warnings.filterwarnings('ignore')

PART 1): Reading the RADAR file¶

In [2]:

filename = "CAPPI-2km-Depth-H&V.nc"

In [3]:

grid = xr.open_dataset(filename)
display(grid)

<xarray.Dataset>
Dimensions:                     (time: 1, x: 500, y: 500, z: 1, nradar: 1)
Coordinates:
  * time                        (time) datetime64[ns] 2022-02-15T21:00:11
  * x                           (x) float64 -2.5e+05 -2.49e+05 ... 2.5e+05
  * y                           (y) float64 -2.5e+05 -2.49e+05 ... 2.5e+05
  * z                           (z) float64 2e+03
Dimensions without coordinates: nradar
Data variables: (12/16)
    origin_latitude             (time) float64 ...
    origin_longitude            (time) float64 ...
    origin_altitude             (time) float64 ...
    projection                  int32 ...
    ProjectionCoordinateSystem  int32 ...
    radar_latitude              (nradar) float64 ...
    ...                          ...
    DBZV                        (time, z, y, x) float32 ...
    UV                          (time, z, y, x) float32 ...
    UH                          (time, z, y, x) float32 ...
    DBZH                        (time, z, y, x) float32 ...
    R_KDP_Corr_h_depth          (time, z, y, x) float32 ...
    R_KDP_Corr_v_depth          (time, z, y, x) float32 ...
Attributes: (12/14)
    Conventions:      Cf/Radial instrument_parameters radar_parameters
    title:            PPIVol
    institution:      EEC
    references:       EEC
    source:           EDGE
    history:          original
    ...               ...
    site_name:        9921GUA
    n_gates_vary:     false
    volume_number:    1
    platform_type:    fixed
    instrument_type:  radar
    primary_axis:     axis_z

PART 2): Reading rain gauges¶

In [4]:

# Load data from text file
obs = np.loadtxt('obs.txt')
obs.shape

Out[4]:

(85,)

In [5]:

#obs

In [6]:

# Load coordinate text file data
obs_coords = np.loadtxt('obs_coords.txt')
obs_coords.shape

Out[6]:

(85, 2)

In [7]:

#obs_coords

PART 3):Preparing the data¶

In [8]:

# select base 2d grid
grid = grid.isel(time=0, z=0)

Radar Horizontal¶

create mask
fill nan with zero
stack xy into radar_obs

In [9]:

nanmask_h = xr.where(np.isnan(grid.R_KDP_Corr_h_depth), True, False)
radar_h = grid.R_KDP_Corr_h_depth.fillna(0)
radar_h_stack = radar_h.stack(radar_obs=("x", "y"))
display(radar_h_stack)

Radar Vertical¶

create mask
fill nan with zero
stack xy into radar_obs

In [10]:

nanmask_v = xr.where(np.isnan(grid.R_KDP_Corr_v_depth), True, False)
radar_v = grid.R_KDP_Corr_v_depth.fillna(0)
radar_v.plot()
radar_v_stack = radar_v.stack(radar_obs=("x", "y"))
display(radar_v_stack)

PART 3): Adjustment of Rain Gauges**¶

In [11]:

# get grid coords
grid_coords = xr.concat([radar_h_stack.x, radar_h_stack.y], "xy").transpose(..., "xy")
display(grid_coords)

AdjustADD with OrdinaryKriging¶

In [12]:

# Estrutura dos seus dados
print("obs:", obs.shape)
print("obs_coords:", obs_coords.shape)
print("radar_h:", radar_h.shape)
print("grid_coords:", grid_coords.shape)
print("radar_v:", radar_v.shape)

obs: (85,)
obs_coords: (85, 2)
radar_h: (500, 500)
grid_coords: (250000, 2)
radar_v: (500, 500)

In [13]:

%%time
# Additive Error Model with OrdinaryKriging Exponential
# class wradlib.ipol.OrdinaryKriging(src, trg, cov='1.0 Exp(10000.)', nnearest=12)
addadjuster_ok_exp = wrl.adjust.AdjustAdd(obs_coords, 
                                      grid_coords.values,
                                      cov = '1.0 Exp(10000.)',
                                      ipclass=wrl.ipol.OrdinaryKriging)
addadjusted_ok_exp_values = addadjuster_ok_exp(obs, radar_h_stack.values)

CPU times: user 20.9 s, sys: 243 ms, total: 21.2 s
Wall time: 21 s

create DataArray, set values, unstack and transpose¶

In [14]:

addadjusted_ok_exp = xr.zeros_like(radar_h_stack)
addadjusted_ok_exp.values = addadjusted_ok_exp_values
addadjusted_ok_exp = addadjusted_ok_exp.unstack().transpose("y", "x")
display(addadjusted_ok_exp)

In [15]:

# open figure
fig = plt.figure(figsize=(10, 3.5))

# "Radar Rainfall: Real"
ax = fig.add_subplot(121)
radar_h.where(~nanmask_h).plot(ax=ax)
ax.set_title("Radar Rainfall: Real")

# "Additive Error Model with Nearest"
ax = fig.add_subplot(122)
addadjusted_ok_exp.where(~nanmask_h).plot(ax=ax)
ax.set_title("Additive Error Model with OK EXPl")

plt.tight_layout()
# salva figura

AdjustADD with ExternalDriftKriging¶

Docstring excerpt:

src_drift : :class:`numpy:numpy.ndarray`
        ndarray of floats, shape (nsrcpoints, )
        values of the external drift at each source point
trg_drift : :class:`numpy:numpy.ndarray`
        ndarray of floats, shape (ntrgpoints, )
        values of the external drift at each target point

Not sure, if the now provided values for src_drift and trg_drift are the correct ones to provide there as drift.

In [16]:

%%time
addadjuster_kde_exp = wrl.adjust.AdjustAdd(obs_coords, 
                                           grid_coords.values,
                                           cov = '1.0 Exp(10000.)',
                                           nnearest=12,
                                           src_drift=obs, # need obs values here
                                           trg_drift=radar_v_stack.values, # need radar v values here
                                           ipclass=wrl.ipol.ExternalDriftKriging)
addadjusted_kde_exp_values = addadjuster_kde_exp(obs, radar_h_stack.values)

CPU times: user 22.8 s, sys: 237 ms, total: 23 s
Wall time: 22.9 s

create DataArray, set values, unstack and transpose¶

In [17]:

addadjusted_kde_exp = xr.zeros_like(radar_h_stack)
addadjusted_kde_exp.values = addadjusted_kde_exp_values
addadjusted_kde_exp = addadjusted_kde_exp.unstack().transpose("y", "x")
display(addadjusted_kde_exp)

In [18]:

# open figure
fig = plt.figure(figsize=(10, 3.5))

# "Radar Rainfall: Real"
ax = fig.add_subplot(121)
radar_h.where(~nanmask_h).plot(ax=ax)
ax.set_title("Radar Rainfall: Real")

# "Additive Error Model with Nearest"
ax = fig.add_subplot(122)
addadjusted_kde_exp.where(~nanmask_h).plot(ax=ax)
ax.set_title("Additive Error Model with OK EXPl")

plt.tight_layout()
# salva figura