Skip to content

fix bug with 100K grid zone because of using representative point not… #485

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
iferencik merged 1 commit into from
Nov 4, 2025
Merged

fix bug with 100K grid zone because of using representative point not… #485

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 22 additions & 6 deletions rapida/components/landuse/search_utils/mgrsconv.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,7 +20,7 @@ def _normalize_lon(lon: float) -> float:

def utm_grid_zone(longitude: float) -> int:
lon = _normalize_lon(float(longitude))
if math.isclose(lon, 180.0, abs_tol=1e-12):
if lon >= 180.0 - 1e-12:
return 60
return int(math.floor((lon + 180.0 - EPS) / 6.0)) + 1

Expand Down Expand Up @@ -119,16 +119,32 @@ def mgrs_100k_tiles_for_bbox(lon_min: float, lat_min: float, lon_max: float, lat
crs_utm = _utm_crs(zone, north)
part_utm = _project_geom(hemi_geom, crs_wgs, crs_utm)

# Build 100k grid cells and label using representative point INSIDE the in-zone part
# Build 100k grid cells and label using geometric center
for cell in generate_utm_100k_grid(part_utm):
inter = cell.intersection(part_utm)
if inter.is_empty:
continue
p = inter.representative_point()
lon_c, lat_c = _wgs84_of_point(p.x, p.y, zone, north)
key = _mgrs_100k_key_for_zone(lat_c, lon_c, zone) # force zone; band+100k from dd2mgrs

# Use the exact geometric center of the 100k cell (always inside the cell)
xmin, ymin, xmax, ymax = cell.bounds
cx = 0.5 * (xmin + xmax)
cy = 0.5 * (ymin + ymax)

# Convert center to WGS84 to get band/100k via your existing helper
lon_c, lat_c = _wgs84_of_point(cx, cy, zone, north)

# Clamp longitude strictly inside THIS zone to avoid seam jitter
lon_min = -180.0 + (zone - 1) * 6.0
lon_max = lon_min + 6.0 - 1e-9
if lon_c < lon_min:
lon_c = lon_min + 1e-9
elif lon_c > lon_max:
lon_c = lon_max

key = _mgrs_100k_key_for_zone(lat_c, lon_c, zone) # still uses latlon2mgrs inside

if key not in out:
out[key] = (box(*cell.bounds), crs_utm) # store full cell polygon (UTM) + CRS
out[key] = (box(*cell.bounds), crs_utm)

return out

Expand Down
51 changes: 28 additions & 23 deletions rapida/components/landuse/search_utils/search.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@ def expand_timerange(start_date: str, end_date: str, days: int = 7) -> tuple[str
return expanded_start.date().isoformat(), expanded_end.date().isoformat()


def process_item(it, mgrs_id, ref_ts, mgrs_poly, mgrs_crs):
def create_candidate(it, mgrs_id, ref_ts, mgrs_poly, mgrs_crs):
props = it.get("properties", {})
item_id = it.get("id", "")
tile_info = get_tileinfo(it)
Expand Down Expand Up @@ -88,6 +88,7 @@ def search( client=None, collection="sentinel-2-l1c",
search_progress = None
cloud_cover = max_cloud_cover
tile_candidates: list[Candidate] = []
item_ids = []
while True:

#this is because the start tiles on every row in every UTM zone do not contain data
Expand All @@ -111,8 +112,12 @@ def search( client=None, collection="sentinel-2-l1c",
"eo:cloud_cover": {"lte": cloud_cover},
}
)

items = [itm.to_dict() for itm in search_result.items()]
if item_ids:
items = [e for e in items if e.get('id') not in item_ids]
item_ids = item_ids.union([e.get('id') for e in items])
else:
item_ids = set([itm.get('id') for itm in items])

logger.debug(
f'Found {len(items)} items in MGRS-{mgrs_id} between {start_date} and {end_date} and {cloud_cover}% cloud cover')
Expand All @@ -125,7 +130,8 @@ def search( client=None, collection="sentinel-2-l1c",
if stop.is_set():break

with ThreadPoolExecutor(max_workers=10) as executor:
futures = [executor.submit(process_item, itm, mgrs_id, mid, mgrs_poly, crs) for itm in items]

futures = [executor.submit(create_candidate, itm, mgrs_id, mid, mgrs_poly, crs) for itm in items]
try:
for future in as_completed(futures):
cand = future.result()
Expand All @@ -134,16 +140,8 @@ def search( client=None, collection="sentinel-2-l1c",
if cand.id in ids:
logger.debug(f'Skipping duplicate tile {cand.id}')
else:
if tile_candidates:
last_cand = tile_candidates[-1]

last_cand_norm = last_cand.tile_data_geometry.intersection(last_cand.mgrs_geometry)
cand_norm = cand.tile_data_geometry.intersection(cand.mgrs_geometry)
similarity = last_cand_norm.intersection(cand_norm).area / last_cand_norm.union(cand_norm).area * 100
if similarity<50:
tile_candidates.append(cand)
else:
tile_candidates.append(cand)
tile_candidates.append(cand)

except KeyboardInterrupt:
for future in futures:
if not future.cancelled():
Expand All @@ -153,17 +151,20 @@ def search( client=None, collection="sentinel-2-l1c",

#sort candidates
scandidates = sorted(tile_candidates, key=lambda c:-c.quality_score)

#merge/union the polygons covering the valid pixels until the whole MGRS 100K tile is covered
union = unary_union([x.tile_data_geometry for x in scandidates]).intersection(scandidates[0].mgrs_geometry)
#union = union.intersection(scandidates[0].mgrs_geometry)
#coverage = int(floor(union.area/mgrs_poly.area*100))
coverage = union.area/mgrs_poly.area*100
union = unary_union([x.tile_data_geometry for x in scandidates])
# normalize by intersecting with ideal MGRS 100k poly
norm_union = union.intersection(mgrs_poly)
#compute coverage
coverage = norm_union.area/mgrs_poly.area*100

# exit here because the whole tile can be covered
if math.isclose(coverage, b=100, abs_tol=1e-2):
tile_candidates = scandidates
logger.debug(f'Found {len(tile_candidates)} suitable candidates in {mgrs_id} between {start_date} and {end_date}')
break
if len(scandidates)==6:break
# no suitable candidates found, enlarge the time interval by one week on both sides and try again
start_date, end_date = expand_timerange(start_date=start_date, end_date=end_date, )
if stop.is_set(): break
Expand Down Expand Up @@ -202,7 +203,6 @@ def search( client=None, collection="sentinel-2-l1c",
cov_poly = c.tile_data_geometry.intersection(mgrs_poly).union(cov_poly)
coverage = cov_poly.area/mgrs_poly.area*100
pruned.append(c)

if math.isclose(coverage, b=100, abs_tol=1e-2):
break
return pruned
Expand Down Expand Up @@ -275,7 +275,6 @@ def fetch_s2_tiles(
description=f'[red]Found {len(cands)} Sentinel2 {image_word} in MGRS-{gid}',
advance=1)
except Exception as e:
print(e)
failed[gid] = e
if progress is not None and total_task is not None:
progress.update(total_task,advance=1)
Expand Down Expand Up @@ -304,6 +303,7 @@ def fetch_s2_tiles(

from rapida.util.setup_logger import setup_logger
from rapida.components.landuse.search_utils.s2item import Sentinel2Item

# logger.setLevel(logging.DEBUG)
logger = setup_logger(level=logging.INFO)
CATALOG_URL = "https://earth-search.aws.element84.com/v1"
Expand All @@ -314,15 +314,20 @@ def fetch_s2_tiles(
UGAKEN_BBOX = 33.280908600000004, -1.154692598710874, 36.59833289999998, 2.650255597059965


bbox = 33, -1, 37, 2



with Progress() as progress:
bbox=UGAKEN_BBOX
bbox=bbox

max_cloud_cover=3
start_date = "2024-03-01"
end_date = "2024-03-30"
max_cloud_cover=2
start_date = "2025-07-01"
end_date = "2025-08-30"
results = fetch_s2_tiles(bbox=bbox,stac_url=CATALOG_URL,
start_date=start_date, end_date=end_date,
max_cloud_cover=max_cloud_cover, progress=progress, prune=True)
max_cloud_cover=max_cloud_cover, progress=progress,filter_for_dev=['36NWF'])


bands = ['B01', 'B02', 'B03', 'B04', 'B05']
Expand Down