Skip to content

Commit

Permalink
Merge
Browse files Browse the repository at this point in the history
  • Loading branch information
@ghiggi
ghiggi committed Feb 5, 2025
2 parents 2536690 + d3a4639 commit 064f05d
Show file tree
Hide file tree
Showing 8 changed files with 1,015 additions and 671 deletions.
6 changes: 6 additions & 0 deletions gpm/accessor/methods.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,12 @@ def crop_around_point(self, lon, lat, distance=None, size=None):

return crop_around_point(self._obj, lon=lon, lat=lat, distance=distance, size=size)

@auto_wrap_docstring
def crop_around_valid_data(self, variable=None):
from gpm.utils.manipulations import crop_around_valid_data

return crop_around_valid_data(self._obj, variable=variable)

@auto_wrap_docstring
def get_crop_slices_by_extent(self, extent):
from gpm.utils.geospatial import get_crop_slices_by_extent
Expand Down
279 changes: 152 additions & 127 deletions gpm/gv/routines.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,7 +36,7 @@
from shapely import Point

import gpm
from gpm.gv.plots import calibration_summary, compare_maps
from gpm.gv.plots import calibration_summary, plot_gdf_map
from gpm.utils.manipulations import (
conversion_factors_degree_to_meter,
)
Expand Down Expand Up @@ -403,87 +403,124 @@ def convert_s_to_ku_band(ds_gr, bright_band_height, z_variable="DBZH"):
return da_ku


def plot_quicklook(
ds_sr,
ds_gr,
min_gr_range,
max_gr_range,
z_min_threshold_gr,
z_min_threshold_sr,
z_variable="DBZH",
cmap="Spectral_r",
):

# Display nearSurface SR reflectivity
da_sr = ds_sr["zFactorFinal"].gpm.slice_range_at_bin(ds_sr["binClutterFreeBottom"])

# Mask by specified sensitivity
mask_gr = ds_gr[z_variable] > z_min_threshold_gr
mask_sr = da_sr > z_min_threshold_sr
def add_radar_info(ax, ds_gr, radar_size):
# - Add radar location
ax.scatter(0, 0, c="black", marker="X", s=radar_size)
ax.scatter(0, 0, c="black", marker="X", s=radar_size)

# Retrieve SR extent
sr_extent = ds_sr.gpm.extent() # noqa
gr_extent = ds_gr.xradar_dev.extent()
ds_gr.xradar_dev.plot_range_distance(
distance=15_000,
ax=ax,
add_background=False,
add_gridlines=False,
add_labels=False,
linestyle="dashed",
edgecolor="black",
)
ds_gr.xradar_dev.plot_range_distance(
distance=100_000,
ax=ax,
add_background=False,
add_gridlines=False,
add_labels=False,
linestyle="dashed",
edgecolor="black",
)
ds_gr.xradar_dev.plot_range_distance(
distance=150_000,
ax=ax,
add_background=False,
add_gridlines=False,
add_labels=False,
linestyle="dashed",
edgecolor="black",
)

# Create figure
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8, 4), dpi=300, subplot_kw={"projection": ccrs.PlateCarree()})

# Plot SR
da_sr.where(mask_sr).gpm.plot_map(
ax=ax1,
cmap=cmap,
vmin=0,
vmax=40,
cbar_kwargs={"label": "SR Surface Reflectivity (dBz)"},
add_colorbar=True,
def plot_quicklook(ds_gr, gdf, sr_z_column, gr_z_column):
# Define Cartopy projection
ccrs_gr_aeqd = ccrs.AzimuthalEquidistant(
central_longitude=ds_gr["longitude"].item(),
central_latitude=ds_gr["latitude"].item(),
)
# Plot GR
ds_gr[z_variable].where(mask_gr).xradar_dev.plot_map(
ax=ax2,
cmap=cmap,
vmin=0,
vmax=40,
cbar_kwargs={"label": "GR Sweep Reflectivity (dBz)", "extend": "both"},
add_colorbar=True,
subplot_kwargs = {}
subplot_kwargs["projection"] = ccrs_gr_aeqd

# Define geographic extent
extent_xy = gdf.total_bounds[[0, 2, 1, 3]]

# Retrieve plot kwargs
plot_kwargs, cbar_kwargs = gpm.get_plot_kwargs("zFactorFinal", user_plot_kwargs={"vmin": 15, "vmax": 45})

# Define figure settings
figsize = (8, 4)
dpi = 300

# Define radar markersize
radar_size = 40

# Create the figure
fig, axes = plt.subplots(1, 3, width_ratios=[1, 1, 1.1], subplot_kw=subplot_kwargs, figsize=figsize, dpi=dpi)

#### Plot SR data
axes[0].coastlines()
_ = plot_gdf_map(
ax=axes[0],
gdf=gdf,
column=sr_z_column,
title="SR Matched",
extent_xy=extent_xy,
# Gridline settings
# grid_linewidth=grid_linewidth,
# grid_color=grid_color,
# Colorbar settings
add_colorbar=False,
# Plot settings
cbar_kwargs=cbar_kwargs,
**plot_kwargs,
)
# - Add the SR swath boundary
da_sr.gpm.plot_swath_lines(ax=ax2)
# - Restrict the extent to the SR overpass
# ax1.set_extent(sr_extent)
# ax2.set_extent(sr_extent)
ax1.set_extent(gr_extent)
ax2.set_extent(gr_extent)
# - Display GR range distances
for ax in [ax1, ax2]:
ds_gr.xradar_dev.plot_range_distance(
distance=min_gr_range,
ax=ax,
add_background=True,
linestyle="dashed",
edgecolor="black",
)
ds_gr.xradar_dev.plot_range_distance(
distance=max_gr_range,
ax=ax,
add_background=True,
linestyle="dashed",
edgecolor="black",
)
ds_gr.xradar_dev.plot_range_distance(
distance=250_000,
ax=ax,
add_background=True,
linestyle="dashed",
edgecolor="black",
add_radar_info(ax=axes[0], ds_gr=ds_gr, radar_size=radar_size)

#### - Plot GR matched data
axes[1].coastlines()
_ = plot_gdf_map(
ax=axes[1],
gdf=gdf,
column=gr_z_column,
title="GR Matched",
extent_xy=extent_xy,
# Gridline settings
# grid_linewidth=grid_linewidth,
# grid_color=grid_color,
# Colorbar settings
add_colorbar=False,
# Plot settings
cbar_kwargs=cbar_kwargs,
**plot_kwargs,
)
add_radar_info(ax=axes[1], ds_gr=ds_gr, radar_size=radar_size)

#### - Plot GR sweep data
axes[2].coastlines()
p = (
ds_gr["DBZH"]
.where(ds_gr["DBZH"] > 0)
.xradar_dev.plot_map(
ax=axes[2],
x="x",
y="y",
add_background=False,
add_gridlines=False,
add_labels=False,
add_colorbar=True,
cbar_kwargs=cbar_kwargs,
**plot_kwargs,
)
# - Add GR location
ax1.scatter(ds_gr["longitude"], ds_gr["latitude"], c="black", marker="X", s=4)
ax2.scatter(ds_gr["longitude"], ds_gr["latitude"], c="black", marker="X", s=4)
# - Set title
ax1.set_title("GPM", fontsize=12, loc="left")
ax2.set_title("Ground Radar", fontsize=12, loc="left")
# - Improve layout and display
plt.tight_layout()
)
p.axes.set_xlim(extent_xy[0:2])
p.axes.set_ylim(extent_xy[2:4])
p.axes.set_title("GR PPI")
add_radar_info(ax=axes[2], ds_gr=ds_gr, radar_size=radar_size)
return fig


Expand Down Expand Up @@ -618,6 +655,7 @@ def volume_matching(
sr_sensitivity_thresholds=None,
download_sr=True,
display_quicklook=True,
display_calibration_summary=False,
quicklook_fpath=None,
):
"""
Expand All @@ -638,11 +676,11 @@ def volume_matching(
beamwidth_gr : float, optional
The GR beamwidth. The default is 1.0.
z_min_threshold_gr : float, optional
The minimum reflectivity threshold (in dBZ) for GR. The default is 0.
Values below this threshold are set to np.nan before matching SR/GR gates.
The minimum reflectivity threshold (in dBZ) for GR. The default is 0 dBZ.
This is used to discard GR gates with ``z_variable_gr`` reflectivity below the threshold.
z_min_threshold_sr : float, optional
The minimum reflectivity threshold (in dBZ) for SR. The default is 10.
Values below this threshold are set to np.nan before matching SR/GR gates.
The minimum reflectivity threshold (in dBZ) for SR. The default is 10 dBZ.
This is used to discard SR footprints with zMeasured reflectivity below the threshold.
min_gr_range : float, optional
The minimum GR range distance (in meters) to select for matching SR/GR gates. The default is 0.
max_gr_range : float, optional
Expand Down Expand Up @@ -732,6 +770,16 @@ def volume_matching(
# - Following Morris and Schwaller 2011 recommendation
ds_gr = ds_gr.where(ds_gr[z_variable_gr] >= z_min_threshold_gr)

#### - Crop GR only to area with data
try:
ds_gr = ds_gr.gpm.crop_around_valid_data(variable=z_variable_gr)
for dim, size in ds_gr.sizes.items():
if size <= 1:
raise ValueError(f"Only a single {dim} has valid data.")
except Exception:
# No valid data for GR
return None

#### - Retrieve GR extent (in WGS84)
extent_gr = ds_gr.xradar_dev.extent(max_distance=None, crs=None)

Expand Down Expand Up @@ -767,30 +815,6 @@ def volume_matching(
# Put SR data into memory
ds_sr = ds_sr.compute()

# Compute SR attenuation correction
ds_sr["zFactorCorrection"] = ds_sr["zFactorFinal"] - ds_sr["zFactorMeasured"]

# Mask SR gates below clutter free region
ds_sr = ds_sr.gpm.mask_below_bin(bins=ds_sr["binClutterFreeBottom"] + 1, strict=False)

####-----------------------------------------------------------------------------.
#### Plot Quicklook
# - TODO: remove because mislading
# --> surface vs radar sweeep scanning high into atmosphere
# if display_quicklook:
# fig = plot_quicklook(
# ds_sr=ds_sr,
# ds_gr=ds_gr,
# min_gr_range=min_gr_range,
# max_gr_range=max_gr_range,
# z_min_threshold_gr=z_min_threshold_gr,
# z_min_threshold_sr=z_min_threshold_sr,
# z_variable=z_variable_gr,
# )
# if quicklook_fpath is not None:
# fig.savefig(quicklook_fpath)
# plt.show()

####-----------------------------------------------------------------------------.
#### Retrieve SR/GR gate resolution, volume and coordinates
#### - Retrieve GR gate resolution
Expand Down Expand Up @@ -844,6 +868,13 @@ def volume_matching(
ds_sr["flagPrecipitationType"] = ds_sr.gpm.retrieve("flagPrecipitationType", method="major_rain_type")
ds_sr["flagHydroClass"] = ds_sr.gpm.retrieve("flagHydroClass")

# Mask SR gates below clutter free region
# - flagHydroClass defines the gates affected by Clutter
ds_sr = ds_sr.gpm.mask_below_bin(bins=ds_sr["binClutterFreeBottom"] + 1, strict=False)

# Compute SR attenuation correction (with masked clutter)
ds_sr["zFactorCorrection"] = ds_sr["zFactorFinal"] - ds_sr["zFactorMeasured"]

####-----------------------------------------------------------------------------.
#### Identify SR gates intersecting GR sweep
# - Define mask of SR footprints within GR range
Expand All @@ -862,7 +893,7 @@ def volume_matching(
####-----------------------------------------------------------------------------.
#### Define SR mask
# Select above minimum SR reflectivity
# mask_sr_minimum_z = ds_sr["zFactorFinal"] > z_sr_min_threshold
mask_sr_minimum_z = ds_sr["zFactorMeasured"] >= z_min_threshold_sr

# Select only 'high quality' data
# - qualityFlag derived from qualityData.
Expand All @@ -882,9 +913,9 @@ def volume_matching(

# Define 3D mask of SR gates matching GR PPI gates
# - TIP: do not mask eccessively here ... but rather at final stage
mask_matched_ppi_3d = mask_ppi & mask_sr_precip & mask_sr_quality
mask_matched_ppi_3d = mask_ppi & mask_sr_precip & mask_sr_quality & mask_sr_minimum_z

# Define 2D mask of SR beams matching the GR PPI
# Define 2D mask of SR rainy beams matching the GR PPI
mask_matched_ppi_2d = mask_matched_ppi_3d.any(dim="range")

####-----------------------------------------------------------------------------.
Expand Down Expand Up @@ -985,7 +1016,7 @@ def volume_matching(
# Initialize Dataset where to add aggregated SR gates
ds_sr_match_ppi = xr.Dataset()

# Mask SR gates not matching the GR PPI
# Mask SR beams not matching the GR PPI and not rainy
ds_sr_ppi = ds_sr[sr_variables].where(mask_matched_ppi_3d)

# Compute aggregation statistics
Expand Down Expand Up @@ -1266,27 +1297,21 @@ def volume_matching(
gdf_match["SR_z_lower_bound"] = gdf_match["SR_z_min"] - gdf_match["SR_vres_mean"] / 2
gdf_match["SR_z_upper_bound"] = gdf_match["SR_z_max"] + gdf_match["SR_vres_mean"] / 2

####-----------------------------------------------------------------------------.
#### Display matching results
if display_quicklook:
# Compare matched volumes
####----------------------------------------------------------------------.
#### Create Quicklook Figure
if display_quicklook or quicklook_fpath:
sr_z_column = f"SR_zFactorFinal_{radar_band}_mean"
gr_z_column = "GR_Z_mean"
fig = compare_maps(
gdf_match,
sr_column=sr_z_column,
gr_column=gr_z_column,
sr_label="SR Reflectivity (dBz)",
gr_label="GR Reflectivity (dBz)",
cmap="Spectral_r",
unified_color_scale=True,
vmin=12,
vmax=40,
)
fig.tight_layout()
plt.show()

# Display raw calibration summary
fig = plot_quicklook(ds_gr=ds_gr, gdf=gdf_match, sr_z_column=sr_z_column, gr_z_column=gr_z_column)
if quicklook_fpath is not None:
fig.savefig(quicklook_fpath)
plt.close()
else:
plt.show()

####----------------------------------------------------------------------.
#### Create Calibration Summary Figure
if display_calibration_summary:
calibration_summary(
df=gdf_match,
gr_z_column=gr_z_column,
Expand Down
Loading

0 comments on commit 064f05d

Please sign in to comment.