Skip to content

Commit

Permalink
Merge pull request #11 from umr-lops/landmask-management
Browse files Browse the repository at this point in the history
keep all tiles , NaN in wind direction variables for tiles with landmask==True
  • Loading branch information
agrouaze authored Feb 15, 2024
2 parents 4a17e91 + 08e56a4 commit 73bbec2
Show file tree
Hide file tree
Showing 4 changed files with 130 additions and 122 deletions.
96 changes: 52 additions & 44 deletions l2awinddirection/generate_L2A_winddir_pdf_product.py
Original file line number Diff line number Diff line change
Expand Up @@ -38,101 +38,109 @@ def generate_wind_distribution_product(tiles, m64rn4, nb_classes=36, shape=(44,
tiles_stacked = tiles.stack(all_tiles=["burst", "tile_line", "tile_sample"])
except:
tiles_stacked = tiles.stack(all_tiles=["tile_line", "tile_sample"])
tiles_stacked_no_nan = tiles_stacked.where(
~np.any(np.isnan(tiles_stacked.sigma0), axis=(0, 1)), drop=True
)
X = tiles_stacked_no_nan.sigma0.transpose("all_tiles", "azimuth", "range").values
# mask0 = ~np.any(
# np.isnan(tiles_stacked.sigma0), axis=(0, 1)
# ) # original mask from Robin, sigma0 never contains NaN contrarily to sigma0_filt
mask = tiles_stacked.land_flag == False # test AG
mask_indices = np.where(mask)
tiles_stacked_no_nan = tiles_stacked.where(mask, drop=True)

if "sigma0_filt" in tiles_stacked:
X = tiles_stacked_no_nan.sigma0_filt.transpose(
"all_tiles", "azimuth", "range"
).values
else:
X = tiles_stacked_no_nan.sigma0.transpose(
"all_tiles", "azimuth", "range"
).values

X_normalized = np.array(
[((x - np.average(x)) / np.std(x)).reshape(shape) for x in X]
) # Normalize data

heading_angle = np.deg2rad(tiles_stacked_no_nan["ground_heading"].values)

heading_angle = np.deg2rad(tiles_stacked["ground_heading"].values)
# Assign new coordinates 'bin_centers' to the dataset
dx = 180 / nb_classes
angles = np.arange(0 + dx / 2, 180 + dx / 2, dx)

predictions = np.ones((tiles_stacked.longitude.size, angles.size)) * np.nan
tiles_stacked_no_nan = tiles_stacked_no_nan.assign_coords(
bin_centers=("bin_centers", angles)
)

# Predict wind direction distributions and add them to the dataset
# predictions = np.squeeze(m64rn4.model.predict(X_normalized)) # original Robin
predictions = m64rn4.model.predict(X_normalized)
tiles_stacked_no_nan = tiles_stacked_no_nan.assign(
if X_normalized.size > 0:
predictions_usable = m64rn4.model.predict(X_normalized)
predictions[mask_indices, :] = predictions_usable

tiles_stacked = tiles_stacked.assign(
wind_direction_distribution=(["all_tiles", "bin_centers"], predictions)
)

# Compute mean wind direction, most likely wind direction and associated standard deviation
mean_wdir = xr.apply_ufunc(
compute_mean_direction,
tiles_stacked_no_nan.wind_direction_distribution,
tiles_stacked_no_nan.bin_centers,
tiles_stacked.wind_direction_distribution,
tiles_stacked.bin_centers,
input_core_dims=[["bin_centers"], ["bin_centers"]],
vectorize=True,
)

mean_wdir = mean_wdir + heading_angle
most_likely_wdir = xr.apply_ufunc(
compute_most_likely_direction,
tiles_stacked_no_nan.wind_direction_distribution,
tiles_stacked_no_nan.bin_centers,
tiles_stacked.wind_direction_distribution,
tiles_stacked.bin_centers,
input_core_dims=[["bin_centers"], ["bin_centers"]],
vectorize=True,
)

most_likely_wdir = most_likely_wdir + heading_angle
std_wdir = xr.apply_ufunc(
compute_standard_deviation,
tiles_stacked_no_nan.wind_direction_distribution,
tiles_stacked_no_nan.bin_centers,
tiles_stacked.wind_direction_distribution,
tiles_stacked.bin_centers,
input_core_dims=[["bin_centers"], ["bin_centers"]],
vectorize=True,
)

# Format dataset
tiles_stacked_no_nan = tiles_stacked_no_nan.assign(
mean_wdir=(["all_tiles"], mean_wdir.data)
)
tiles_stacked_no_nan["mean_wdir"].attrs[
tiles_stacked = tiles_stacked.assign(mean_wdir=(["all_tiles"], mean_wdir.data))
tiles_stacked["mean_wdir"].attrs[
"definition"
] = "180° ambiguous mean wind direction."
tiles_stacked_no_nan["mean_wdir"].attrs["long_name"] = "Mean wind direction"
tiles_stacked_no_nan["mean_wdir"].attrs[
tiles_stacked["mean_wdir"].attrs["long_name"] = "Mean wind direction"
tiles_stacked["mean_wdir"].attrs[
"convention"
] = "Clockwise, relative to geographic North."
tiles_stacked_no_nan["mean_wdir"].attrs["units"] = "degree"
tiles_stacked_no_nan["mean_wdir"].attrs["vmin"] = 0
tiles_stacked_no_nan["mean_wdir"].attrs["vmax"] = 180
tiles_stacked["mean_wdir"].attrs["units"] = "degree"
tiles_stacked["mean_wdir"].attrs["vmin"] = 0
tiles_stacked["mean_wdir"].attrs["vmax"] = 180

tiles_stacked_no_nan = tiles_stacked_no_nan.assign(
tiles_stacked = tiles_stacked.assign(
most_likely_wdir=(["all_tiles"], most_likely_wdir.data)
)
tiles_stacked_no_nan["most_likely_wdir"].attrs[
tiles_stacked["most_likely_wdir"].attrs[
"definition"
] = "180° ambiguous most likely wind direction, defined with a bin precision of 2.5°."
tiles_stacked_no_nan["most_likely_wdir"].attrs[
"long_name"
] = "Most likely wind direction"
tiles_stacked_no_nan["most_likely_wdir"].attrs[
tiles_stacked["most_likely_wdir"].attrs["long_name"] = "Most likely wind direction"
tiles_stacked["most_likely_wdir"].attrs[
"convention"
] = "Clockwise, relative to geographic North."
tiles_stacked_no_nan["most_likely_wdir"].attrs["units"] = "degree"
tiles_stacked_no_nan["most_likely_wdir"].attrs["vmin"] = 0
tiles_stacked_no_nan["most_likely_wdir"].attrs["vmax"] = 180
tiles_stacked["most_likely_wdir"].attrs["units"] = "degree"
tiles_stacked["most_likely_wdir"].attrs["vmin"] = 0
tiles_stacked["most_likely_wdir"].attrs["vmax"] = 180

tiles_stacked_no_nan = tiles_stacked_no_nan.assign(
std_wdir=(["all_tiles"], std_wdir.data)
)
tiles_stacked_no_nan["std_wdir"].attrs[
tiles_stacked = tiles_stacked.assign(std_wdir=(["all_tiles"], std_wdir.data))
tiles_stacked["std_wdir"].attrs[
"definition"
] = "Standard deviation associated with wind direction."
tiles_stacked_no_nan["std_wdir"].attrs[
tiles_stacked["std_wdir"].attrs[
"long_name"
] = "Standard deviation of the wind direction"
tiles_stacked_no_nan["std_wdir"].attrs["units"] = "degree"
tiles_stacked_no_nan["std_wdir"].attrs["vmin"] = 0
tiles_stacked_no_nan["std_wdir"].attrs["vmax"] = 180
tiles_stacked["std_wdir"].attrs["units"] = "degree"
tiles_stacked["std_wdir"].attrs["vmin"] = 0
tiles_stacked["std_wdir"].attrs["vmax"] = 180

return tiles_stacked_no_nan.unstack()
return tiles_stacked.unstack()


def compute_mean_direction(weights, angles):
Expand Down
130 changes: 60 additions & 70 deletions l2awinddirection/generate_L2A_winddir_regression_product.py
Original file line number Diff line number Diff line change
Expand Up @@ -2,20 +2,18 @@
inspired from predict_wind_direction.py (project_rmarquart)
January 2024
"""
import pdb
import numpy as np
import scipy as sp
import xarray as xr
from scipy import stats
import glob
import os, sys
from tqdm import tqdm
from l2awinddirection.M64RN4 import M64RN4_regression


# --------------------- #
# ----- Functions ----- #
# --------------------- #

def generate_wind_product(tiles, model_regs, shape = (44, 44, 1)):

def generate_wind_product(tiles, model_regs, shape=(44, 44, 1)):
"""
Generate wind direction from previously extracted tiles and store them in an xarray.Dataset.
Parameters:
Expand All @@ -25,35 +23,63 @@ def generate_wind_product(tiles, model_regs, shape = (44, 44, 1)):
Returns:
(xarray.Dataset): dataset containing the tiles with associated wind direction.
"""
tiles_stacked = tiles.stack(all_tiles = ['burst', 'tile_line','tile_sample'])
tiles_stacked_no_nan = tiles_stacked.where(~np.any(np.isnan(tiles_stacked.sigma0), axis = (0, 1)), drop = True)
X = tiles_stacked_no_nan.sigma0.transpose('all_tiles', 'azimuth', 'range').values
X_normalized = np.array([((x - np.average(x))/np.std(x)).reshape(shape) for x in X]) # Normalize data

heading_angle = np.deg2rad(tiles_stacked_no_nan['ground_heading'].values)

#Predict wind directions using all models and obtain mean and std values
predictions = launch_prediction(X_normalized, model_regs)
mean_prediction = sp.stats.circmean(predictions, np.pi, axis = 1)
std_prediction = sp.stats.circstd(predictions, np.pi, axis = 1)

try:
tiles_stacked = tiles.stack(all_tiles=["burst", "tile_line", "tile_sample"])
except:
tiles_stacked = tiles.stack(all_tiles=["tile_line", "tile_sample"])

mask = tiles_stacked.land_flag == False # test AG
mask_indices = np.where(mask)
tiles_stacked_no_nan = tiles_stacked.where(mask, drop=True)

if "sigma0_filt" in tiles_stacked:
X = tiles_stacked_no_nan.sigma0_filt.transpose(
"all_tiles", "azimuth", "range"
).values
else:
X = tiles_stacked_no_nan.sigma0.transpose(
"all_tiles", "azimuth", "range"
).values

X_normalized = np.array(
[((x - np.average(x)) / np.std(x)).reshape(shape) for x in X]
) # Normalize data

heading_angle = np.deg2rad(tiles_stacked["ground_heading"].values)
predictions = np.ones((tiles_stacked.longitude.size,len(model_regs))) * np.NaN
# Predict wind directions using all models and obtain mean and std values
if X_normalized.size > 0:
predictions_usable = launch_prediction(X_normalized, model_regs)
predictions[mask_indices,:] = predictions_usable

mean_prediction = sp.stats.circmean(predictions, np.pi, axis=1)
std_prediction = sp.stats.circstd(predictions, np.pi, axis=1)

predicted_wdir = mean_prediction + heading_angle

# Format dataset
tiles_stacked_no_nan = tiles_stacked_no_nan.assign(wind_direction = (['all_tiles'], np.rad2deg(predicted_wdir)%180))
tiles_stacked_no_nan['wind_direction'].attrs['definition'] = '180° ambiguous wind direction, clockwise, relative to geographic North.'
tiles_stacked_no_nan['wind_direction'].attrs['units'] = 'degree'
tiles_stacked_no_nan['wind_direction'].attrs['vmin'] = 0
tiles_stacked_no_nan['wind_direction'].attrs['vmax'] = 180

tiles_stacked_no_nan = tiles_stacked_no_nan.assign(wind_std = (['all_tiles'], np.rad2deg(std_prediction)%180))
tiles_stacked_no_nan['wind_std'].attrs['definition'] = 'Standard deviation associated with wind direction. Calculated by computing the spread accross the '
tiles_stacked_no_nan['wind_std'].attrs['units'] = 'degree'
tiles_stacked_no_nan['wind_std'].attrs['vmin'] = 0
tiles_stacked_no_nan['wind_std'].attrs['vmax'] = 180

return tiles_stacked_no_nan.unstack()

tiles_stacked = tiles_stacked.assign(
wind_direction=(["all_tiles"], np.rad2deg(predicted_wdir) % 180)
)
tiles_stacked["wind_direction"].attrs[
"definition"
] = "180° ambiguous wind direction, clockwise, relative to geographic North."
tiles_stacked["wind_direction"].attrs["units"] = "degree"
tiles_stacked["wind_direction"].attrs["vmin"] = 0
tiles_stacked["wind_direction"].attrs["vmax"] = 180

tiles_stacked = tiles_stacked.assign(
wind_std=(["all_tiles"], np.rad2deg(std_prediction) % 180)
)
tiles_stacked["wind_std"].attrs[
"definition"
] = "Standard deviation associated with wind direction. Calculated by computing the spread accross the "
tiles_stacked["wind_std"].attrs["units"] = "degree"
tiles_stacked["wind_std"].attrs["vmin"] = 0
tiles_stacked["wind_std"].attrs["vmax"] = 180

return tiles_stacked.unstack()


def launch_prediction(X_normalized, model_regs):
"""
Expand All @@ -68,41 +94,5 @@ def launch_prediction(X_normalized, model_regs):

for i, m64rn4 in enumerate(model_regs):
predictions[:, i] = np.squeeze(m64rn4.model.predict(X_normalized))

return predictions


if __name__ == "__main__":

#Parameters to instantiante the models
input_shape = (44, 44, 1)
data_augmentation = True
learning_rate = 1e-3

model_regs = []

#Modify path_best_models depending of the acquisition mode (iw, ew, wv)
path_best_models = glob.glob('.../analysis/s1_data_analysis/project_rmarquar/wsat/trained_models/iw/*.hdf5')
for path in path_best_models:

m64rn4_reg = M64RN4_regression(input_shape, data_augmentation)
m64rn4_reg.create_and_compile(learning_rate)

m64rn4_reg.model.load_weights(path)

model_regs.append(m64rn4_reg)


#Get listing of files from which to generate the wind direction
files = glob.glob('.../analysis/s1_data_analysis/project_rmarquar/slc_wind_product/1.4k/*/tiles*1.4k.nc')
print("Number of files to process: %d" % len(files))

for file in tqdm(files):
tiles = xr.open_dataset(file)
if len(tiles.variables.keys()) == 0:
continue
res = generate_wind_product(tiles, model_regs)
res.to_netcdf(file.replace('.nc', '_wind.nc'))
# if you want to remove the file containing the tiles used for inference (data kept in final product)
del tiles
os.remove(file)
return predictions
Original file line number Diff line number Diff line change
Expand Up @@ -37,6 +37,12 @@
# logger.setLevel(logging.FATAL)

import l2awinddirection

try:
print(l2awinddirection.__version__)
print(l2awinddirection.__file__)
except:
pass
from l2awinddirection.generate_L2A_winddir_pdf_product import (
generate_wind_distribution_product,
)
Expand Down Expand Up @@ -200,7 +206,10 @@ def main():
if not os.path.exists(safefile) and args.skipmoves is True:
logging.info(" step 1: move %s -> %s", l2awindirtilessafe, safefile)
shutil.copytree(l2awindirtilessafe, safefile)
files = glob.glob(os.path.join(safefile, "*.nc"))
polarization_usable = "vv"
files = sorted(
glob.glob(os.path.join(safefile, "*" + polarization_usable + "*.nc"))
)
logging.info("Number of files to process: %d" % len(files))
if not os.path.exists(args.outputdir):
logging.info("mkdir %s", args.outputdir)
Expand All @@ -225,9 +234,6 @@ def main():
outputfile = outputfile.replace(workspace, args.outputdir)
if not os.path.exists(os.path.dirname(outputfile)):
os.makedirs(os.path.dirname(outputfile))
res.attrs["winddirection_L2A_processor"] = "l2awinddirection"
res.attrs["winddirection_L2A_processor_version"] = l2awinddirection.__version__
res.to_netcdf(outputfile)
elif args.mode == "regression":

outputfile = file.replace(
Expand All @@ -239,9 +245,13 @@ def main():
os.makedirs(os.path.dirname(outputfile))

res = generate_wind_product(tiles, model_m64rn4)
res.attrs["winddirection_L2A_processor"] = "l2awinddirection"
res.attrs["winddirection_L2A_processor_version"] = l2awinddirection.__version__
res.to_netcdf(outputfile)
res.attrs["winddirection_L2A_processor"] = "l2awinddirection"
res.attrs["winddirection_L2A_processor_version"] = l2awinddirection.__version__
res.attrs["used_channel"] = res.coords["pol"].values
res.attrs["winddirection_L2A_processor_mode"] = args.mode
if 'pol' in res.coords:
res = res.drop('pol')
res.to_netcdf(outputfile)

# if you want to remove the file containing the tiles used for inference (data kept in final product)
del tiles
Expand Down
2 changes: 1 addition & 1 deletion pyproject.toml
Original file line number Diff line number Diff line change
Expand Up @@ -47,4 +47,4 @@ default_section = "THIRDPARTY"
known_first_party = "l2awinddirection"

[project.scripts]
L2a-wind-dir-processor = "l2awinddirection.l2awinddirection:main"
L2a-wind-dir-processor = "l2awinddirection.mainl2awinddirection:main"

0 comments on commit 73bbec2

Please sign in to comment.