Skip to content
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.

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

Add region selection from file #981

Merged
merged 8 commits into from
Sep 18, 2023
Merged
Changes from 1 commit
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
Prev Previous commit
Next Next commit
remove coordinate opt
  • Loading branch information
acordonez committed Sep 12, 2023
commit d7d1c4cfae05e95feb75cead2fbabbc99676dba1
51 changes: 16 additions & 35 deletions pcmdi_metrics/io/region_from_file.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,50 +7,31 @@
import xarray as xr
import xcdat

def region_from_file(data,feature,coords=None,rgn_path=None,attr=None):
# Return data masked from coordinate list or shapefile.
# Masks a single region
def region_from_file(data,rgn_path,attr,feature):
# Return data masked from a feature in input file.
# Arguments:
# data: xcdat dataset
# feature: str, name of region
# coords: list, coordinates
# rgn_path: str, path to file
# attr: str, attribute name

lon = data["lon"].data
lat = data["lat"].data

# Option 1: Region is defined by coord pairs
if coords is not None:
print("Using coordinates to select region.")
try:
names=[feature]
regions = regionmask.Regions([np.array(coords)],names=names)
mask = regions.mask(lon, lat)
val=0
except Exception as e:
print("Error in extracting region by provided coordinates:")
raise e

# Option 2: region is defined by shapefile
elif rgn_path is not None:
print("Reading region from file.")
try:
regions_df = gpd.read_file(rgn_path)
if attr is not None:
regions = regionmask.from_geopandas(regions_df,names=attr)
else:
print("Attribute name not provided.")
regions = regionmask.from_geopandas(regions_df)
mask = regions.mask(lon, lat)
# Can't match mask by name, rather index of name
val = list(regions_file[attr]).index(feature)
except Exception as e:
print("Error in creating region subset from file:")
raise e

else:
raise RuntimeError("Error in region selection: Region coordinates or region file must be provided.")
print("Reading region from file.")
try:
regions_df = gpd.read_file(rgn_path)
if attr is not None:
regions = regionmask.from_geopandas(regions_df,names=attr)
else:
print("Attribute name not provided.")
regions = regionmask.from_geopandas(regions_df)
mask = regions.mask(lon, lat)
# Can't match mask by name, rather index of name
val = list(regions_file[attr]).index(feature)
except Exception as e:
print("Error in creating region subset from file:")
raise e

try:
masked_data = data.where(mask == val)
Expand Down