Skip to content

Commit

Permalink
Merge pull request #12 from agrouaze/fixdeps
Browse files Browse the repository at this point in the history
changes with production scripts WV L1C (logs+args)
  • Loading branch information
agrouaze authored May 23, 2023
2 parents 2c7c5b4 + da68c03 commit 878a3a3
Show file tree
Hide file tree
Showing 7 changed files with 297 additions and 50 deletions.
3 changes: 2 additions & 1 deletion docs/conf.py
Original file line number Diff line number Diff line change
Expand Up @@ -61,7 +61,8 @@
"dask": ("https://docs.dask.org/en/latest", None),
"xarray": ("https://docs.xarray.dev/en/latest/", None),
"rasterio": ("https://rasterio.readthedocs.io/en/latest/", None),
"datatree": ("https://xarray-datatree.readthedocs.io/en/latest/", None)
"datatree": ("https://xarray-datatree.readthedocs.io/en/latest/", None),
'geoviews': ('https://geoviews.org/index.html', None)
}


Expand Down
3 changes: 2 additions & 1 deletion slcl1butils/config.yaml
Original file line number Diff line number Diff line change
Expand Up @@ -2,4 +2,5 @@ ecmwf0.1_pattern: './%Y/%j/ECMWF_FORECAST_0100_%Y%m%d%H%M_10U_10V.nc'
ww3_pattern: './%Y/FIELD_NC/LOPS_WW3-GLOB-30M_%Y%m.nc'
l1c_iw_version: '0.1'
iw_outputdir: none
data_dir: '../assests/'
data_dir: '../assests/'
wv_outputdir: none
84 changes: 49 additions & 35 deletions slcl1butils/get_polygons_from_l1b.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,6 +4,7 @@
from shapely import wkt
import xarray as xr
import logging
from xsarslc.tools import xndindex

polygons_varnames = ['swath', 'tiles', 'bursts']
def get_swath_tiles_polygons_from_l1bgroup(l1b_ds, swath_only=False, ik=0,burst_type='intra', **kwargs):
Expand Down Expand Up @@ -41,7 +42,19 @@ def get_swath_tiles_polygons_from_l1bgroup(l1b_ds, swath_only=False, ik=0,burst_
return polygons, coordinates, variables
logging.debug('l1b_ds : %s',l1b_ds)
# Tiles polynom & variables
for iburst in l1b_ds['burst'].values:

gege = xndindex(l1b_ds['sigma0'].sizes) # whatever tthere is or not tile_line and tile_sample and burst, it loops on it
for uu in gege:
iburst = uu['burst']
if 'tile_line' in uu:
itile_line = uu['tile_line']
else:
itile_line = np.nan
if 'tile_sample' in uu:
itile_sample = uu['tile_sample']
else:
itile_sample = np.nan
#for iburst in l1b_ds['burst'].values:
if burst_type == 'intra':
if 'burst_corner_longitude' in l1b_ds:
burst_lon_corners = l1b_ds['burst_corner_longitude'].sel(burst=iburst).values.flatten().tolist()
Expand All @@ -54,40 +67,41 @@ def get_swath_tiles_polygons_from_l1bgroup(l1b_ds, swath_only=False, ik=0,burst_
#pts_burst = [pts_burst[0], pts_burst[1], pts_burst[3], pts_burst[2]]
#poly_burst = geometry.Polygon(pts_burst)
poly_bursts.append(poly_burst)
for itile_sample in l1b_ds['tile_sample'].values:
for itile_line in l1b_ds['tile_line'].values:
# Find lon/lat tile corners
lon_corners = l1b_ds['corner_longitude'].sel(burst=iburst, tile_sample=itile_sample,
tile_line=itile_line).values.flatten().tolist()
lat_corners = l1b_ds['corner_latitude'].sel(burst=iburst, tile_sample=itile_sample,
tile_line=itile_line).values.flatten().tolist()
# Check on the lon/lat corners validity
# print(np.sum((~np.isfinite(lon_corners))))
if np.sum((~np.isfinite(lon_corners))) == 0:
# Define the tile polygons
pts = [geometry.Point(lon_corners[cpt], lat_corners[cpt]) for cpt, _ in enumerate(lon_corners)]
pts = [pts[0], pts[1], pts[3], pts[2],pts[0]]
logging.debug('pts: %s',pts)
logging.debug('one pt : %s %s',pts[0],type(pts[0]))
order = [0,1,3,2,0]
poly_tile = geometry.Polygon(np.stack([np.array(lon_corners)[order],np.array(lat_corners)[order]]).T)
#poly_tile = geometry.Polygon(pts)
poly_tiles.append(poly_tile)
# Coordinates
ibursts.append(iburst)
itile_samples.append(itile_sample)
itile_lines.append(itile_line)
if variable_names:
if (variable_names[0] is not None):
for variable_name in variable_names:
if (variable_name == 'macs_Im'):
variables[variable_name].append(float(
l1b_ds[variable_name].sel(burst=iburst, tile_sample=itile_sample,
tile_line=itile_line).values[ik]))
else:
variables[variable_name].append(float(
l1b_ds[variable_name].sel(burst=iburst, tile_sample=itile_sample,
tile_line=itile_line).values))
# for itile_sample in l1b_ds['tile_sample'].values:
# for itile_line in l1b_ds['tile_line'].values:
# Find lon/lat tile corners
# lon_corners = l1b_ds['corner_longitude'].sel(burst=iburst, tile_sample=itile_sample,
# tile_line=itile_line).values.flatten().tolist()
# lat_corners = l1b_ds['corner_latitude'].sel(burst=iburst, tile_sample=itile_sample,
# tile_line=itile_line).values.flatten().tolist()
lon_corners = l1b_ds['corner_longitude'].sel(uu).values.flatten().tolist()
lat_corners = l1b_ds['corner_latitude'].sel(uu).values.flatten().tolist()

# Check on the lon/lat corners validity
# print(np.sum((~np.isfinite(lon_corners))))
if np.sum((~np.isfinite(lon_corners))) == 0:
# Define the tile polygons
pts = [geometry.Point(lon_corners[cpt], lat_corners[cpt]) for cpt, _ in enumerate(lon_corners)]
pts = [pts[0], pts[1], pts[3], pts[2],pts[0]]
logging.debug('pts: %s',pts)
logging.debug('one pt : %s %s',pts[0],type(pts[0]))
order = [0,1,3,2,0]
poly_tile = geometry.Polygon(np.stack([np.array(lon_corners)[order],np.array(lat_corners)[order]]).T)
#poly_tile = geometry.Polygon(pts)
poly_tiles.append(poly_tile)
# Coordinates
ibursts.append(iburst)
itile_samples.append(itile_sample)
itile_lines.append(itile_line)
if variable_names:
if (variable_names[0] is not None):
for variable_name in variable_names:
if (variable_name == 'macs_Im'):
variables[variable_name].append(float(
l1b_ds[variable_name].sel(uu).values[ik]))
else:
variables[variable_name].append(float(
l1b_ds[variable_name].sel(uu).values))

polygons['tiles'] = poly_tiles
polygons['bursts'] = poly_bursts
Expand Down
166 changes: 166 additions & 0 deletions slcl1butils/plotting/wallpaper.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,166 @@
import xarray as xr
import numpy as np
import matplotlib.pyplot as plt
import datatree
import os
import logging
from matplotlib import colors as mcolors
cmap = mcolors.LinearSegmentedColormap.from_list("", ["white","violet","mediumpurple","cyan","springgreen","yellow","red"])
PuOr = mcolors.LinearSegmentedColormap.from_list("", ["darkgoldenrod","white","purple"])

def plot_wallpaper(l1b_path, png_path='wallpaper_xspectra.png'):
"""
Args:
l1b_path: str
png_path: str
Returns:
"""
root, ext = os.path.splitext(png_path)
real_path = root+'_real'+ext
imag_path = root+'_imag'+ext
xsplot, orbit_pass = make_wallpaper(l1b_path)
if orbit_pass.upper()=='ASCENDING':
xsplot = xsplot.stack(trange=['tile_sample','range'])
xsplot = xsplot.assign_coords({'burst':np.flip(xsplot['burst'].data), 'tile_line':np.flip(xsplot['tile_line'].data)}).sortby('burst', 'tile_line')
xsplot = xsplot.stack(tazimuth=['burst','tile_line','azimuth'])
elif orbit_pass.upper()=='DESCENDING':
xsplot = xsplot.assign_coords({'tile_sample':np.flip(xsplot['tile_sample'].data)}).sortby('tile_sample')
xsplot = xsplot.stack(trange=['tile_sample','range'])
xsplot = xsplot.stack(tazimuth=['burst','tile_line','azimuth'])
else:
raise ValueError('Unknown orbit pass: {}'.format(orbit_pass))

plt.figure(figsize=(320,180))
plt.imshow(xsplot['xspectra_real'].transpose('tazimuth','trange','rgb').data)
plt.savefig(real_path, dpi='figure',bbox_inches='tight')
plt.close()
logging.info('real path png: %s',real_path)
plt.figure(figsize=(320,180))
plt.imshow(xsplot['xspectra_imag'].transpose('tazimuth','trange','rgb').data)
plt.savefig(imag_path, dpi='figure',bbox_inches='tight')
plt.close()
logging.info('imaginary path png: %s', imag_path)


def make_wallpaper(l1b_path):
"""
Args:
l1b_path: str
Returns:
"""
from xsarslc.processing.xspectra import symmetrize_xspectrum
from xsarslc.tools import xndindex
from matplotlib.backends.backend_agg import FigureCanvasAgg as FigureCanvas

l1b = datatree.open_datatree(l1b_path)
xs2tau = l1b['intraburst']['xspectra_2tau_Re']+1j*l1b['intraburst']['xspectra_2tau_Im']
xs2tau = xs2tau.assign_coords({'k_rg':xs2tau['k_rg'].mean(dim=set(xs2tau['k_rg'].dims)-set(['freq_sample']), keep_attrs=True)}).swap_dims({'freq_sample':'k_rg', 'freq_line':'k_az'})
xs2tau = symmetrize_xspectrum(xs2tau).squeeze(dim='2tau')

xsplot = list()
for mytile in xndindex({'burst':xs2tau.sizes['burst'], 'tile_line':xs2tau.sizes['tile_line'], 'tile_sample':xs2tau.sizes['tile_sample']}):

heading = np.radians(l1b['intraburst']['ground_heading'][mytile].data)
incidence = np.round(l1b['intraburst']['incidence'][mytile].item(),2)
tau = np.round(l1b['intraburst']['tau'][mytile].item(),3)
try:
cutoff = int(l1b['intraburst'].ds[mytile]['azimuth_cutoff'].data)
except:
cutoff = np.NaN
lon = np.round(l1b['intraburst'].ds[mytile]['longitude'].item(), 2)
lat = np.round(l1b['intraburst'].ds[mytile]['latitude'].item(), 2)
keast = xs2tau['k_rg']*np.cos(heading)+xs2tau['k_az']*np.sin(heading)
knorth = xs2tau['k_az']*np.cos(heading)-xs2tau['k_rg']*np.sin(heading)
keast.attrs.update({'long_name':'wavenumber in East direction', 'units':'rad/m'})
knorth.attrs.update({'long_name':'wavenumber in North direction', 'units':'rad/m'})
xs2tau = xs2tau.assign_coords({'k_east':keast,'k_north':knorth})
heading=np.arctan2(np.sin(heading), np.cos(heading))
range_rotation = -np.degrees(heading) if np.abs(np.degrees(heading))<=90 else -np.degrees(heading)+180
azimuth_rotation = -np.degrees(heading)+90 if np.degrees(heading)>=0 else -np.degrees(heading)-90

figreal = plt.figure(figsize=(10,10), tight_layout=True)
xs2tau[mytile].real.plot(cmap=cmap, vmin=0, x='k_east', y='k_north')
ax = plt.gca()
for r in [400,200,100,50]:
circle = plt.Circle((0, 0), 2*np.pi/r, color='k', fill=False, linestyle='--', linewidth=0.5)
ax.add_patch(circle)
plt.text(-np.sqrt(2)*np.pi/r-0.002,np.sqrt(2)*np.pi/r+0.002,'{} m'.format(r), rotation=45.,horizontalalignment='center',verticalalignment='center')
for a in [-60,-30,30,60]:
plt.plot([-0.2*np.cos(np.radians(a)), 0.2*np.cos(np.radians(a))],[-0.2*np.sin(np.radians(a)), 0.2*np.sin(np.radians(a))], color='k', linestyle='--', linewidth=0.5)
plt.vlines(0,-0.2,0.2, color='k', linestyle='--', linewidth=0.5)
plt.hlines(0,-0.2,0.2, color='k', linestyle='--', linewidth=0.5)
xp = 0.14*np.cos(heading)
yp = 0.14*np.sin(heading)
plt.plot([-xp,xp], [yp,-yp], color='r') # range line
plt.plot([yp,-yp], [xp,-xp], color='r') # azimuth line
plt.plot(np.array([-xp,xp])+2*np.pi/cutoff*np.sin(heading), np.array([yp,-yp])+2*np.pi/cutoff*np.cos(heading), color='k', linestyle='--') # cutoff upper line
plt.plot(np.array([-xp,xp])-2*np.pi/cutoff*np.sin(heading), np.array([yp,-yp])-2*np.pi/cutoff*np.cos(heading), color='k', linestyle='--') # cutoff lower line
plt.axis('scaled')
plt.xlim([-0.14,0.14])
plt.ylim([-0.14,0.14])

plt.text(0.9*xp,-0.9*yp+0.006,'range', color='r', rotation=range_rotation, fontsize=15,horizontalalignment='center',verticalalignment='center') # range name
plt.text(0.85*yp+0.006,0.85*xp,'azimuth', color='r', rotation=azimuth_rotation, fontsize=15,horizontalalignment='center',verticalalignment='center') # azimuth name
plt.text(-0.135,-0.135,'cuto-off : {} m'.format(cutoff))
plt.text(-0.135,0.13,'incidence : {} deg'.format(incidence))
plt.text(0.065,0.13,'longitude : {} deg'.format(lon))
plt.text(0.065,0.12,'latitude : {} deg'.format(lat))
plt.text(-0.135,0.12,'tau : {} s'.format(tau))
plt.title('X-spectrum (real)')

canvas = FigureCanvas(figreal)
width, height = figreal.get_size_inches() * figreal.get_dpi()
canvas.draw()
imagereal = np.frombuffer(canvas.tostring_rgb(), dtype='uint8').reshape(int(height), int(width), 3)[146:854,102:810,:]
plt.close()

# -----------------------------------------

figimag = plt.figure(figsize=(10,10), tight_layout=True)
xs2tau[mytile].imag.plot(cmap=PuOr, x='k_east', y='k_north')
ax = plt.gca()
for r in [400,200,100,50]:
circle = plt.Circle((0, 0), 2*np.pi/r, color='k', fill=False, linestyle='--', linewidth=0.5)
ax.add_patch(circle)
plt.text(-np.sqrt(2)*np.pi/r-0.002,np.sqrt(2)*np.pi/r+0.002,'{} m'.format(r), rotation=45.,horizontalalignment='center',verticalalignment='center')
for a in [-60,-30,30,60]:
plt.plot([-0.2*np.cos(np.radians(a)), 0.2*np.cos(np.radians(a))],[-0.2*np.sin(np.radians(a)), 0.2*np.sin(np.radians(a))], color='k', linestyle='--', linewidth=0.5)
plt.vlines(0,-0.2,0.2, color='k', linestyle='--', linewidth=0.5)
plt.hlines(0,-0.2,0.2, color='k', linestyle='--', linewidth=0.5)

plt.plot([-xp,xp], [yp,-yp], color='r') # range line
plt.plot([yp,-yp], [xp,-xp], color='r') # azimuth line
plt.plot(np.array([-xp,xp])+2*np.pi/cutoff*np.sin(heading), np.array([yp,-yp])+2*np.pi/cutoff*np.cos(heading), color='k', linestyle='--') # cutoff upper line
plt.plot(np.array([-xp,xp])-2*np.pi/cutoff*np.sin(heading), np.array([yp,-yp])-2*np.pi/cutoff*np.cos(heading), color='k', linestyle='--') # cutoff lower line

plt.axis('scaled')
plt.xlim([-0.14,0.14])
plt.ylim([-0.14,0.14])
# plt.grid()
plt.text(0.9*xp,-0.9*yp+0.006,'range', color='r', rotation=range_rotation, fontsize=15,horizontalalignment='center',verticalalignment='center') # range name
plt.text(0.85*yp+0.006,0.85*xp,'azimuth', color='r', rotation=azimuth_rotation, fontsize=15,horizontalalignment='center',verticalalignment='center') # azimuth name
plt.text(-0.135,-0.135,'cuto-off : {} m'.format(cutoff))
plt.text(-0.135,0.13,'incidence : {} deg'.format(incidence))
plt.text(-0.135,0.12,'tau : {} s'.format(tau))
plt.text(0.065,0.13,'longitude : {} deg'.format(lon))
plt.text(0.065,0.12,'latitude : {} deg'.format(lat))
plt.title('X-spectrum (imaginary)')

canvas = FigureCanvas(figimag)
width, height = figimag.get_size_inches() * figimag.get_dpi()
canvas.draw()
imageimag = np.frombuffer(canvas.tostring_rgb(), dtype='uint8').reshape(int(height), int(width), 3)[146:854,102:810,:]
plt.close()

xsreal = xr.DataArray(imagereal, dims=('azimuth','range','rgb'), name='xspectra_real').assign_coords(mytile)
xsimag = xr.DataArray(imageimag, dims=('azimuth','range','rgb'), name='xspectra_imag').assign_coords(mytile)
xs = xr.merge([xsreal, xsimag])
xsplot.append(xs)
xsplot = xr.combine_by_coords([xs.expand_dims(['burst', 'tile_line','tile_sample']) for xs in xsplot])
return xsplot, l1b['intraburst'].attrs['orbit_pass']
4 changes: 2 additions & 2 deletions slcl1butils/scripts/do_WV_L1C_SAFE_from_L1B_SAFE.pbs
Original file line number Diff line number Diff line change
Expand Up @@ -8,7 +8,7 @@ echo start
# inputs
#$1 str l1bsafe
#$2 str version of the run: 1.5 for instance

#$3 outputdir
myvariable=$(whoami)
. /appli/anaconda/latest/etc/profile.d/conda.sh
#conda activate /home/datawork-cersat-public/cache/project/mpc-sentinel1/workspace/mamba/envs/xsar_oct22
Expand All @@ -18,4 +18,4 @@ which python

exe=do_WV_L1C_SAFE_from_L1B_SAFE

$exe --overwrite --l1bsafe $1 --version $2
$exe --overwrite --l1bsafe $1 --version $2 --outputdir $3
Loading

0 comments on commit 878a3a3

Please sign in to comment.