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.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Create regionmask.py
  • Loading branch information
acordonez committed Sep 8, 2023
commit da7059a49b53fafc41914e784f7314aa8addaac3
53 changes: 53 additions & 0 deletions pcmdi_metrics/io/regionmask.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,53 @@
import geopandas as gpd
import numpy as np
import os
import pandas as pd
import regionmask
import sys
import xarray as xr
import xcdat

def mask_region(data,name,coords=None,shp_path=None,column=None):
# Return data masked from coordinate list or shapefile.
# Masks a single region

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

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

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

else:
raise RuntimeError("Error in masking: Region coordinates or shapefile must be provided.")

try:
masked_data = data.where(mask == val)
except Exception as e:
print("Error: Masking failed.")
raise e

return masked_data