Skip to content

Commit

Permalink
update to add dimension if only not exists and add new way to transla…
Browse files Browse the repository at this point in the history
…te shapefile to 360
  • Loading branch information
sliu008 committed Nov 20, 2023
1 parent ea6668f commit fcd9898
Show file tree
Hide file tree
Showing 2 changed files with 35 additions and 4 deletions.
3 changes: 2 additions & 1 deletion podaac/subsetter/group_handling.py
Original file line number Diff line number Diff line change
Expand Up @@ -123,7 +123,8 @@ def recombine_grouped_datasets(datasets: List[xr.Dataset], output_file: str, sta
for dim_name in list(dataset.dims.keys()):
new_dim_name = dim_name.split(GROUP_DELIM)[-1]
dim_group = _get_nested_group(base_dataset, dim_name)
dim_group.createDimension(new_dim_name, dataset.dims[dim_name])
if new_dim_name not in dim_group.dimensions:
dim_group.createDimension(new_dim_name, dataset.dims[dim_name])

# Rename variables
_rename_variables(dataset, base_dataset, start_date, time_vars)
Expand Down
36 changes: 33 additions & 3 deletions podaac/subsetter/subset.py
Original file line number Diff line number Diff line change
Expand Up @@ -962,9 +962,39 @@ def subset_with_shapefile(dataset: xr.Dataset,
def convert_180_to_360(lon, lat):
return tuple(map(lambda value: value + 360 if value < 0 else value, lon)), lat

geometries = [transform(convert_180_to_360, geometry) for geometry in
shapefile_df.geometry]
shapefile_df.geometry = geometries
def translate_longitude(geometry):
def translate_point(point):
# Translate the point's x-coordinate (longitude) by adding 360
return Point((point.x + 360) % 360, point.y)

def translate_polygon(polygon):
def translate_coordinates(coords):
first_coord = coords[0]
if len(coords[0]) == 2:
return [((x+360)%360, y) for x,y in coords]
elif len(coords[0]) == 3:
return [((x+360)%360, y, z) for x,y,z in coords]

exterior = translate_coordinates(polygon.exterior.coords)

interiors = [
translate_coordinates(ring.coords)
for ring in polygon.interiors
]

return Polygon(exterior, interiors)

if isinstance(geometry, (Point, Polygon)):
return translate_point(geometry) if isinstance(geometry, Point) else translate_polygon(geometry)
elif isinstance(geometry, MultiPolygon):
# Translate each polygon in the MultiPolygon
translated_polygons = [translate_longitude(subgeometry) for subgeometry in geometry]
return MultiPolygon(translated_polygons)
else:
# Handle other geometry types as needed
return geometry

shapefile_df.geometry = shapefile_df['geometry'].apply(translate_longitude)

# Mask and scale shapefile
def scale(lon, lat):
Expand Down

0 comments on commit fcd9898

Please sign in to comment.