Skip to content

Commit 3f077c2

Browse files
konstntokaskonstntokas
committed
backup
1 parent 42c9869 commit 3f077c2

2 files changed

Lines changed: 242 additions & 39 deletions

File tree

‎examples/sen1_analysismode/debug_resampling_sen3_new.py‎

Lines changed: 62 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -2,10 +2,70 @@
22
from rectify_new import rectify_dataset
33
import matplotlib.pyplot as plt
44
import numpy as np
5+
import dask
6+
7+
import contextily as ctx
8+
import geopandas as gpd
9+
from shapely.geometry import Point
10+
11+
12+
def plot_points(xs, ys, ax):
13+
# --- Build GeoDataFrame
14+
gdf = gpd.GeoDataFrame(
15+
geometry=[Point(lon, lat) for lon, lat in zip(xs, ys)],
16+
crs="EPSG:4326", # lat/lon
17+
)
18+
19+
# --- Reproject to Web Mercator (EPSG:3857) for tiles
20+
gdf = gdf.to_crs(epsg=3857)
21+
22+
# --- Plot
23+
24+
# plot the scatter
25+
gdf.plot(ax=ax, color="red", marker="o", label="Your points")
26+
27+
# add base map
28+
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, zoom=7)
29+
30+
ax.set_title("Pixels with zero area")
31+
# ax.set_axis_off()
32+
33+
34+
dask.config.set(scheduler="synchronous")
535

636
source_ds = xr.open_zarr("./data/S3-OLCI-L2A.zarr.zip", consolidated=False)
7-
# print(np.where(np.diff(source_ds.lon.compute(), axis=1) == 0))
8-
# print(np.where(np.diff(source_ds.lat.compute(), axis=0) == 0))
37+
# jx, ix = np.where(np.diff(source_ds.lon.compute(), axis=1) == 0)
38+
# jy, iy = np.where(np.diff(source_ds.lat.compute(), axis=0) == 0)
39+
# j_all = np.concatenate((jx, jy))
40+
# i_all = np.concatenate((ix, iy))
941
target_ds = rectify_dataset(source_ds, interp_methods="bilinear")
1042
target_ds.rtoa_8.plot(vmin=0.0, vmax=0.3)
1143
plt.show()
44+
45+
# fi, ax = plt.subplots(1, 2, figsize=(16, 10))
46+
# target_ds.rtoa_8.plot(ax=ax[0], vmin=0.0, vmax=0.3)
47+
# plot_points(
48+
# source_ds.lon.values[j_all, i_all], source_ds.lat.values[j_all, i_all], ax[1]
49+
# )
50+
# plt.tight_layout()
51+
# plt.show()
52+
53+
# url = "https://objects.eodc.eu/e05ab01a9d56408d82ac32d69a5aae2a:202601-s03slslst-eu/29/products/cpm_v262/S3B_SL_2_LST____20260129T195739_20260129T200039_20260129T223030_0179_116_128_0720_ESA_O_NR_004.zarr"
54+
#
55+
# dt = xr.open_datatree(url, engine="zarr", chunks={})
56+
# source_ds = dt.measurements.to_dataset()
57+
# jx, ix = np.where(np.diff(source_ds.longitude.compute(), axis=1) == 0)
58+
# jy, iy = np.where(np.diff(source_ds.latitude.compute(), axis=0) == 0)
59+
# j_all = np.concatenate((jx, jy))
60+
# i_all = np.concatenate((ix, iy))
61+
# target_ds = rectify_dataset(source_ds, interp_methods="nearest")
62+
#
63+
# fi, ax = plt.subplots(1, 2, figsize=(16, 10))
64+
# target_ds.lst.plot(ax=ax[0])
65+
# plot_points(
66+
# source_ds.longitude.values[j_all, i_all],
67+
# source_ds.latitude.values[j_all, i_all],
68+
# ax[1],
69+
# )
70+
# plt.tight_layout()
71+
# plt.show()

‎examples/sen1_analysismode/rectify_new.py‎

Lines changed: 180 additions & 37 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,15 @@
5252
normalize_grid_mapping,
5353
)
5454

55+
import matplotlib.pyplot as plt
56+
import contextily as ctx
57+
import geopandas as gpd
58+
from shapely.geometry import Point
59+
60+
61+
FX_MIN = FY_MIN = -0.001
62+
FXFY_MAX = 1.002
63+
5564

5665
def rectify_dataset(
5766
source_ds: xr.Dataset,
@@ -753,24 +762,34 @@ def _rectify_block(
753762
target_y_valid = target_y_coords[valid]
754763
x_diff = x_valid - target_x_valid
755764
y_diff = y_valid - target_y_valid
756-
ix_start = np.where(x_diff > 0, ix_valid, ix_valid - 1)
765+
ix_start = np.where(x_diff < 0, ix_valid, ix_valid - 1)
757766
iy_start = np.where(y_diff > 0, iy_valid, iy_valid - 1)
758-
valid_edge = (
759-
(ix_start >= 0)
760-
& (ix_start <= x_coords.shape[1] - 2)
761-
& (iy_start >= 0)
762-
& (iy_start <= y_coords.shape[0] - 2)
767+
x_coords = np.pad(
768+
x_coords,
769+
pad_width=((0, 2), (0, 2)),
770+
mode="constant",
771+
constant_values=np.nan,
772+
)
773+
y_coords = np.pad(
774+
y_coords,
775+
pad_width=((0, 2), (0, 2)),
776+
mode="constant",
777+
constant_values=np.nan,
763778
)
764-
ix_start = ix_start[valid_edge]
765-
iy_start = iy_start[valid_edge]
766-
target_x_valid = target_x_valid[valid_edge]
767-
target_y_valid = target_y_valid[valid_edge]
768779
x_corners = np.column_stack(
769780
(
770781
x_coords[iy_start, ix_start],
771782
x_coords[iy_start, ix_start + 1],
772783
x_coords[iy_start + 1, ix_start],
773784
x_coords[iy_start + 1, ix_start + 1],
785+
x_coords[iy_start, ix_start + 1],
786+
x_coords[iy_start, ix_start + 2],
787+
x_coords[iy_start + 1, ix_start + 1],
788+
x_coords[iy_start + 1, ix_start + 2],
789+
x_coords[iy_start + 1, ix_start],
790+
x_coords[iy_start + 1, ix_start + 1],
791+
x_coords[iy_start + 2, ix_start],
792+
x_coords[iy_start + 2, ix_start + 1],
774793
)
775794
)
776795
y_corners = np.column_stack(
@@ -779,6 +798,14 @@ def _rectify_block(
779798
y_coords[iy_start, ix_start + 1],
780799
y_coords[iy_start + 1, ix_start],
781800
y_coords[iy_start + 1, ix_start + 1],
801+
y_coords[iy_start, ix_start + 1],
802+
y_coords[iy_start, ix_start + 2],
803+
y_coords[iy_start + 1, ix_start + 1],
804+
y_coords[iy_start + 1, ix_start + 2],
805+
y_coords[iy_start, ix_start + 1],
806+
y_coords[iy_start, ix_start + 2],
807+
y_coords[iy_start + 1, ix_start + 1],
808+
y_coords[iy_start + 1, ix_start + 2],
782809
)
783810
)
784811
fx, fy, valid_fracs = bilinear_fractions(
@@ -793,8 +820,8 @@ def _rectify_block(
793820
value_u0 = value_00 + fx * (value_01 - value_00)
794821
value_u1 = value_10 + fx * (value_11 - value_10)
795822
valid0, valid1 = np.where(valid)
796-
valid0 = valid0[valid_edge][valid_fracs]
797-
valid1 = valid1[valid_edge][valid_fracs]
823+
valid0 = valid0[valid_fracs]
824+
valid1 = valid1[valid_fracs]
798825
data_rectified[:, valid0, valid1] = value_u0 + fy * (value_u1 - value_u0)
799826
else:
800827
raise NotImplementedError(
@@ -813,9 +840,9 @@ def bilinear_fractions(
813840
"""Compute bilinear fractions (fx, fy) for targets inside quads.
814841
815842
Args:
816-
x_corners : shape (n, 4) x-coordinates of quad corners,
843+
x_corners : shape (n, 12) x-coordinates of quad corners,
817844
ordered as [p00, p01, p10, p11]
818-
y_corners : shape (n, 4) y-coordinates of quad corners, same order
845+
y_corners : shape (n, 12) y-coordinates of quad corners, same order
819846
x_target : shape (n,) target x coordinates
820847
y_target : shape (n,) target y coordinates
821848
@@ -824,45 +851,161 @@ def bilinear_fractions(
824851
fy : shape (n,) fraction along y-direction
825852
valid_frac: shape (n,) boolean for filtering
826853
"""
827-
assert x_corners.shape[1] == 4
828-
assert y_corners.shape[1] == 4
854+
assert x_corners.shape[1] == 12
855+
assert y_corners.shape[1] == 12
856+
857+
run_idx = np.arange(y_target.size)
858+
859+
# upper left triangle
860+
fx00, fy00, valid_frac00 = _calc_fx_fy(
861+
x_corners[:, :4], y_corners[:, :4], x_target, y_target, run_idx
862+
)
863+
run_idx = run_idx[~np.isin(run_idx, valid_frac00)]
864+
865+
# lower right triangle
866+
fx03, fy03, valid_frac03 = _calc_fx_fy(
867+
x_corners[run_idx, :4],
868+
y_corners[run_idx, :4],
869+
x_target[run_idx],
870+
y_target[run_idx],
871+
run_idx,
872+
p3_corner=True,
873+
)
874+
run_idx = run_idx[~np.isin(run_idx, valid_frac03)]
875+
876+
# upper left triangle x_shift
877+
fx10, fy10, valid_frac10 = _calc_fx_fy(
878+
x_corners[run_idx, 4:8],
879+
y_corners[run_idx, 4:8],
880+
x_target[run_idx],
881+
y_target[run_idx],
882+
run_idx,
883+
)
884+
run_idx = run_idx[~np.isin(run_idx, valid_frac10)]
885+
886+
# lower right triangle x_shift
887+
fx13, fy13, valid_frac13 = _calc_fx_fy(
888+
x_corners[run_idx, 4:8],
889+
y_corners[run_idx, 4:8],
890+
x_target[run_idx],
891+
y_target[run_idx],
892+
run_idx,
893+
p3_corner=True,
894+
)
895+
run_idx = run_idx[~np.isin(run_idx, valid_frac13)]
896+
897+
# upper left triangle y_shift
898+
fx20, fy20, valid_frac20 = _calc_fx_fy(
899+
x_corners[run_idx, 8:12],
900+
y_corners[run_idx, 8:12],
901+
x_target[run_idx],
902+
y_target[run_idx],
903+
run_idx,
904+
)
905+
run_idx = run_idx[~np.isin(run_idx, valid_frac20)]
906+
907+
# lower right triangle x_shift
908+
fx23, fy23, valid_frac23 = _calc_fx_fy(
909+
x_corners[run_idx, 8:12],
910+
y_corners[run_idx, 8:12],
911+
x_target[run_idx],
912+
y_target[run_idx],
913+
run_idx,
914+
p3_corner=True,
915+
)
829916

830-
valid_frac = np.arange(y_target.size)
917+
valid_frac = np.concatenate(
918+
(
919+
valid_frac00,
920+
valid_frac03,
921+
valid_frac10,
922+
valid_frac13,
923+
valid_frac20,
924+
valid_frac23,
925+
)
926+
)
927+
fx = np.concatenate((fx00, fx03, fx10, fx13, fx20, fx23))
928+
fy = np.concatenate((fy00, fy03, fy10, fy13, fy20, fy23))
929+
930+
return fx, fy, valid_frac
931+
932+
933+
def _calc_fx_fy(x_corners, y_corners, x_target, y_target, run_idx, p3_corner=False):
831934

832935
valid_corners = np.all(~np.isnan(x_corners), axis=1)
833936
x_corners = x_corners[valid_corners]
834937
y_corners = y_corners[valid_corners]
835938
x_target = x_target[valid_corners]
836939
y_target = y_target[valid_corners]
837-
valid_frac = valid_frac[valid_corners]
940+
run_idx = run_idx[valid_corners]
941+
942+
if p3_corner:
943+
dx1 = x_corners[:, 2] - x_corners[:, 3]
944+
dx2 = x_corners[:, 1] - x_corners[:, 3]
945+
dy1 = y_corners[:, 2] - y_corners[:, 3]
946+
dy2 = y_corners[:, 1] - y_corners[:, 3]
947+
x0 = x_corners[:, 3]
948+
y0 = y_corners[:, 3]
949+
else:
950+
dx1 = x_corners[:, 1] - x_corners[:, 0]
951+
dx2 = x_corners[:, 2] - x_corners[:, 0]
952+
dy1 = y_corners[:, 1] - y_corners[:, 0]
953+
dy2 = y_corners[:, 2] - y_corners[:, 0]
954+
x0 = x_corners[:, 0]
955+
y0 = y_corners[:, 0]
838956

839-
dx1 = x_corners[:, 1] - x_corners[:, 0]
840-
dx2 = x_corners[:, 2] - x_corners[:, 0]
841-
dy1 = y_corners[:, 1] - y_corners[:, 0]
842-
dy2 = y_corners[:, 2] - y_corners[:, 0]
843957
det = dx1 * dy2 - dx2 * dy1
844958
valid_det = det != 0
845-
846-
# dx1 = x_corners[:, 2] - x_corners[:, 3]
847-
# dx2 = x_corners[:, 1] - x_corners[:, 3]
848-
# dy1 = y_corners[:, 2] - y_corners[:, 3]
849-
# dy2 = y_corners[:, 1] - y_corners[:, 3]
850-
# det_b = dx1 * dy2 - dx2 * dy1
851-
# valid_det = det != 0
852-
853-
x_corners = x_corners[valid_det]
854-
y_corners = y_corners[valid_det]
959+
x0 = x0[valid_det]
960+
y0 = y0[valid_det]
855961
x_target = x_target[valid_det]
856962
y_target = y_target[valid_det]
857963
dx1 = dx1[valid_det]
858964
dx2 = dx2[valid_det]
859965
dy1 = dy1[valid_det]
860966
dy2 = dy2[valid_det]
861967
det = det[valid_det]
862-
valid_frac = valid_frac[valid_det]
968+
run_idx = run_idx[valid_det]
863969

864-
# vectorized solve: fx = ( (tx-x0)*dy2 - (ty-y0)*dx2 ) / det
865-
fx = ((x_target - x_corners[:, 0]) * dy2 - (y_target - y_corners[:, 0]) * dx2) / det
866-
fy = (dx1 * (y_target - y_corners[:, 0]) - dy1 * (x_target - x_corners[:, 0])) / det
970+
fx = ((x_target - x0) * dy2 - (y_target - y0) * dx2) / det
971+
fy = (dx1 * (y_target - y0) - dy1 * (x_target - x0)) / det
972+
valid_corner = (fx >= FX_MIN) & (fy >= FY_MIN) & ((fx + fy) <= FXFY_MAX)
867973

868-
return fx, fy, valid_frac
974+
return fx[valid_corner], fy[valid_corner], run_idx[valid_corner]
975+
976+
977+
def plot_points(xs, ys):
978+
# --- Build GeoDataFrame
979+
gdf = gpd.GeoDataFrame(
980+
geometry=[Point(lon, lat) for lon, lat in zip(xs, ys)],
981+
crs="EPSG:4326", # lat/lon
982+
)
983+
984+
# --- Reproject to Web Mercator (EPSG:3857) for tiles
985+
gdf = gdf.to_crs(epsg=3857)
986+
987+
# --- Calculate reasonable extent around points + margin
988+
bbox = (
989+
gpd.GeoSeries([Point(12, 55), Point(20, 65)], crs="EPSG:4326")
990+
.to_crs(3857)
991+
.total_bounds
992+
)
993+
994+
# --- Plot
995+
fig, ax = plt.subplots(figsize=(8, 10))
996+
997+
# plot the scatter
998+
gdf.plot(ax=ax, color="red", marker="o", label="Your points")
999+
1000+
# set extent (optional override for Southern Sweden)
1001+
ax.set_xlim(bbox[0], bbox[2])
1002+
ax.set_ylim(bbox[1], bbox[3])
1003+
1004+
# add base map
1005+
ctx.add_basemap(ax, source=ctx.providers.OpenStreetMap.Mapnik, zoom=7)
1006+
1007+
ax.set_title("Southern Sweden Scatter with Base Map")
1008+
ax.set_axis_off()
1009+
1010+
plt.legend()
1011+
plt.show()

0 commit comments

Comments
 (0)