From 3d5d159d7fb820f8f696091ccf6fc817fe5ed231 Mon Sep 17 00:00:00 2001 From: xhong Date: Fri, 5 Sep 2025 12:15:34 +0800 Subject: [PATCH] =?UTF-8?q?fix(HLS=5FSuPER):=20=E7=A1=AE=E4=BF=9D=E5=A4=9A?= =?UTF-8?q?=E8=BE=B9=E5=BD=A2=E5=A4=96=E7=8E=AF=E4=B8=BA=E9=80=86=E6=97=B6?= =?UTF-8?q?=E9=92=88=E6=96=B9=E5=90=91=E5=B9=B6=E7=AE=80=E5=8C=96=E7=BD=91?= =?UTF-8?q?=E7=BB=9C=E6=9F=A5=E8=AF=A2=E6=97=B6=E6=8F=90=E4=BA=A4=E7=9A=84?= =?UTF-8?q?=E5=9D=90=E6=A0=87.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- HLS_SuPER/HLS_SuPER.py | 34 +++++++++++++++++++++------------- 1 file changed, 21 insertions(+), 13 deletions(-) diff --git a/HLS_SuPER/HLS_SuPER.py b/HLS_SuPER/HLS_SuPER.py index b3453ab..826b3cc 100644 --- a/HLS_SuPER/HLS_SuPER.py +++ b/HLS_SuPER/HLS_SuPER.py @@ -5,7 +5,7 @@ HLS Subsetting, Processing, and Exporting Reformatted Data Prep Script Authors: Cole Krehbiel, Mahsa Jami, and Erik Bolch Contact: lpdaac@usgs.gov Editor: Hong Xie -Last Updated: 2025-01-15 +Last Updated: 2025-09-05 =============================================================================== """ @@ -24,6 +24,7 @@ import json import earthaccess from shapely.geometry import box +from shapely.geometry.polygon import orient import geopandas as gpd from datetime import datetime as dt import dask.distributed @@ -171,7 +172,15 @@ def parse_arguments(): return parser.parse_args() - +def ensure_ccw(geom): + """ + Ensure the exterior ring of the polygon is counterclockwise. + """ + if geom.exterior.is_ccw: + return geom # Already counterclockwise + else: + return orient(geom, sign=1.0) # Make it counterclockwise + def format_roi(roi): """ Determines if submitted ROI is a file or bbox coordinates. @@ -192,22 +201,20 @@ def format_roi(roi): logging.info( "Multiple polygons detected. Creating single geometry of external coordinates." ) - single_geometry = roi.unary_union.convex_hull + single_geometry = roi.union_all().convex_hull roi = gpd.GeoDataFrame(geometry=[single_geometry], crs=roi.crs) # Check if ROI is in Geographic CRS, if not, convert to it - if roi.crs.is_geographic: - # List Vertices in correct order for search - # (Add) 使用外包矩形坐标作为检索使用的坐标 - minx, miny, maxx, maxy = roi.total_bounds - bounding_box = box(minx, miny, maxx, maxy) - vertices_list = list(bounding_box.exterior.coords) - - else: - roi_geographic = roi.to_crs("EPSG:4326") + if not roi.crs.is_geographic: + roi = roi.to_crs("EPSG:4326") logging.info( "Note: ROI submitted is being converted to Geographic CRS (EPSG:4326)" ) - vertices_list = list(roi_geographic.geometry[0].exterior.coords) + roi['geometry'] = roi['geometry'].apply(ensure_ccw) + # List Vertices in correct order for search + # (Add) 使用外包矩形坐标简化提交检索时使用的坐标 + minx, miny, maxx, maxy = roi.total_bounds + bounding_box = box(minx, miny, maxx, maxy) + vertices_list = list(bounding_box.exterior.coords) except (FileNotFoundError, ValueError): sys.exit( f"The GeoJSON/shapefile is either not valid or could not be found.\nPlease double check the name and provide the absolute path to the file or make sure that it is located in {os.getcwd()}" @@ -219,6 +226,7 @@ def format_roi(roi): # Convert bbox to a geodataframe for clipping roi = gpd.GeoDataFrame(geometry=[box(*bbox)], crs="EPSG:4326") + roi['geometry'] = roi['geometry'].apply(ensure_ccw) vertices_list = list(roi.geometry[0].exterior.coords)