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)