Commit 98e742bd authored by CEVAER's avatar CEVAER
Browse files

Fixed var orientation, do not need to_north_ref anymore

parent 183297ef
......@@ -29,8 +29,7 @@ def save_netcdf(sat_file, path, center_status, percent_outside, percent_inside_i
ds = xr.open_dataset(sat_file_sw)
ds_s = ds.squeeze()
ds_polar = to_polar(ds_s, lon_center, lat_center)
ds_polar_north = to_north_ref(ds_polar)
ds_polar_north = to_polar(ds_s, lon_center, lat_center)
eye_coords = []
for c_i in range(0, len(lons_eye_shape)):
......
......@@ -39,25 +39,16 @@ def custom_meshgrid(theta, rad, data_shape):
return coord_theta, coord_rad
def lon_lat_orientation(polar_ds):
polar_ds.lon.values[:] = np.flip(polar_ds.lon.values, axis=0)
polar_ds.lat.values[:] = np.flip(polar_ds.lat.values, axis=0)
def var_orientation(var, idx_azim, teta_shape):
var[:] = np.flip(var, axis=0)
t_shp = polar_ds.theta.shape
# Get the index in theta array corresponding to the 90° angle. (as theta[idx_90] =~ 90)
idx_90 = int(t_shp[0] / 4)
# Getting azimuth to north for center pixel
azim_heading = polar_ds.heading_angle.sel(theta=0, rad=0)
# Getting the closest index to the azim_heading azimuth
idx_azim = np.searchsorted(polar_ds.theta.values, azim_heading, side="left")
idx_90 = int(teta_shape[0] / 4)
# Roll the lon and lat arrays of 90° + azim_angle to match the correct orientation.
polar_ds.lon.values[:] = np.roll(polar_ds.lon.values, idx_90 + idx_azim, axis=0)
polar_ds.lat.values[:] = np.roll(polar_ds.lat.values, idx_90 + idx_azim, axis=0)
var[:] = np.roll(var, idx_90 + idx_azim, axis=0)
return polar_ds
return var
def to_polar(ds, lon_center, lat_center):
......@@ -75,6 +66,12 @@ def to_polar(ds, lon_center, lat_center):
rads_meter = rads * pix_dist
#print(rads_meter.min(), rads_meter.max())
# Getting azimuth to north for center pixel
azim_heading = ds.heading_angle.isel(x=yloc, y=xloc).item()
# Getting the closest index to the azim_heading azimuth
idx_azim = np.searchsorted(theta, azim_heading, side="left")
coord_theta, coord_rad = np.meshgrid(theta, rads_meter, indexing="ij")
geod = Geod(ellps='WGS84')
......@@ -103,6 +100,7 @@ def to_polar(ds, lon_center, lat_center):
nan_ar[:] = np.nan
# Keeping only valid indices
nan_ar[valid_indices] = var_polar[valid_indices]
nan_ar = var_orientation(nan_ar, idx_azim, theta.shape)
vars_polar[var_name] = (["theta", "rad"], nan_ar)
vars_polar["lon"] = (["theta", "rad"], lons_polar)
......@@ -113,7 +111,7 @@ def to_polar(ds, lon_center, lat_center):
"rad": (["rad"], rads_meter)
})
polar_ds = lon_lat_orientation(polar_ds)
#polar_ds = lon_lat_orientation(polar_ds)
return polar_ds
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment