Skip to content

Commit

Permalink
[pre-commit.ci] auto fixes from pre-commit.com hooks
Browse files Browse the repository at this point in the history
for more information, see https://pre-commit.ci
  • Loading branch information
pre-commit-ci[bot] committed Feb 5, 2024
1 parent a6d90d3 commit 1cd21d3
Show file tree
Hide file tree
Showing 63 changed files with 8,206 additions and 8,011 deletions.
2 changes: 1 addition & 1 deletion carbonplan_forest_risks/__init__.py
Original file line number Diff line number Diff line change
Expand Up @@ -5,5 +5,5 @@
try:
version = get_distribution(__name__).version
except DistributionNotFound: # pragma: no cover
version = '0.0.0' # pragma: no cover
version = "0.0.0" # pragma: no cover
__version__ = version
12 changes: 6 additions & 6 deletions carbonplan_forest_risks/collect.py
Original file line number Diff line number Diff line change
Expand Up @@ -4,11 +4,11 @@

def fire(y, src, inds=None):
da = xr.Dataset()
da['x'] = src.x
da['y'] = src.y
da['time'] = src.time
da['lat'] = (['y', 'x'], src['lat'])
da['lon'] = (['y', 'x'], src['lon'])
da["x"] = src.x
da["y"] = src.y
da["time"] = src.time
da["lat"] = (["y", "x"], src["lat"])
da["lon"] = (["y", "x"], src["lon"])
shape = (len(src.time), len(src.y), len(src.x))

if inds is None:
Expand All @@ -17,6 +17,6 @@ def fire(y, src, inds=None):
yfull = np.zeros(shape).flatten()
yfull[inds] = y
yfull = yfull.reshape(shape)
da['prediction'] = (['time', 'y', 'x'], yfull)
da["prediction"] = (["time", "y", "x"], yfull)

return da
25 changes: 15 additions & 10 deletions carbonplan_forest_risks/fit/growth.py
Original file line number Diff line number Diff line change
Expand Up @@ -14,13 +14,13 @@ def logistic(x, f, p):
)


def growth(x, y, f, noise='gamma', init=None):
def growth(x, y, f, noise="gamma", init=None):
def loglik(x, y, f, p):
a, b, c, w0, w1, scale = p
_mu = logistic(x, f, [a, b, c, w0, w1])
if noise == 'gamma':
if noise == "gamma":
return -np.sum(gamma.logpdf(y, np.maximum(_mu / scale, 1e-12), scale=scale))
if noise == 'normal':
if noise == "normal":
return -np.sum(norm.logpdf(y, loc=_mu, scale=scale))

fx = lambda p: loglik(x, y, f, p)
Expand All @@ -33,14 +33,16 @@ def loglik(x, y, f, p):
bounds = Bounds(lb, ub)
if init is None:
init = [np.nanmean(y), 0.1, 10, 0, 0, np.nanstd(y)]
options_trust = {'maxiter': 10000}
options_trust = {"maxiter": 10000}

with warnings.catch_warnings():
warnings.filterwarnings('ignore', category=UserWarning)
result = minimize(fx, init, bounds=bounds, method='trust-constr', options=options_trust)
warnings.filterwarnings("ignore", category=UserWarning)
result = minimize(
fx, init, bounds=bounds, method="trust-constr", options=options_trust
)

if result.success is False:
print('optimization failed')
print("optimization failed")

return GrowthModel(result, noise, x, y, f)

Expand All @@ -63,12 +65,15 @@ def r2(self, x, f, y):

def predict(self, x, f, percentile=None):
if percentile is not None:
f = [np.nanpercentile(f[0], percentile[0]), np.nanpercentile(f[1], percentile[1])]
f = [
np.nanpercentile(f[0], percentile[0]),
np.nanpercentile(f[1], percentile[1]),
]
return logistic(x, f, self.p)

def sample(self, x, f):
mu = logistic(x, f, self.p)
if self.noise == 'gamma':
if self.noise == "gamma":
return np.random.gamma(mu / self.scale, scale=self.scale)
if self.noise == 'normal':
if self.noise == "normal":
return np.random.normal(mu, scale=self.scale)
13 changes: 10 additions & 3 deletions carbonplan_forest_risks/fit/hurdle.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,11 +9,16 @@ def hurdle(x, y, log=True, max_iter=1000):
x, y = remove_nans(x, y)
n_obs = len(x)

clf = LogisticRegression(fit_intercept=True, penalty='none', max_iter=max_iter)
clf = LogisticRegression(fit_intercept=True, penalty="none", max_iter=max_iter)

if log:
reg = TweedieRegressor(
fit_intercept=True, power=0, link='log', alpha=0, tol=1e-8, max_iter=max_iter
fit_intercept=True,
power=0,
link="log",
alpha=0,
tol=1e-8,
max_iter=max_iter,
)
else:
reg = LinearRegression(fit_intercept=True)
Expand All @@ -31,7 +36,9 @@ def __init__(self, clf, reg, n_obs, log=None, x=None, y=None):
self.log = log
self.n_obs = n_obs
if x is not None and y is not None:
self.train_r2 = np.corrcoef(self.predict_linear(x)[y > 0], y[y > 0])[0, 1] ** 2
self.train_r2 = (
np.corrcoef(self.predict_linear(x)[y > 0], y[y > 0])[0, 1] ** 2
)
self.train_roc = roc_auc_score(y > 0, self.predict_prob(x))

def __repr__(self):
Expand Down
12 changes: 6 additions & 6 deletions carbonplan_forest_risks/fit/interp.py
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,11 @@
def score(y_true, y_pred):
m = ~np.isnan(y_pred)
if len(y_pred) != m.sum():
print(f'found {len(m) - m.sum()} nans in y_pred')
print(f"found {len(m) - m.sum()} nans in y_pred")
return r2_score(y_true[m], y_pred[m])


def interp(df, mask, var='biomass', spacing=4000):
def interp(df, mask, var="biomass", spacing=4000):
"""
Grid a set of lat/lon points to a grid defined by mask
Expand All @@ -38,7 +38,7 @@ def interp(df, mask, var='biomass', spacing=4000):

# extract the projection and grid info
region = [mask.x.data[0], mask.x.data[-1], mask.y.data[-1], mask.y.data[0]]
projection = pyproj.Proj(mask.attrs['crs'])
projection = pyproj.Proj(mask.attrs["crs"])

coordinates = (df.lon.values, df.lat.values)

Expand All @@ -54,8 +54,8 @@ def interp(df, mask, var='biomass', spacing=4000):
# fit the gridder
chain = vd.Chain(
[
('mean', vd.BlockReduce(np.mean, spacing=spacing * 0.25, region=region)),
('nearest', vd.ScipyGridder(method='linear')),
("mean", vd.BlockReduce(np.mean, spacing=spacing * 0.25, region=region)),
("nearest", vd.ScipyGridder(method="linear")),
]
)

Expand All @@ -64,7 +64,7 @@ def interp(df, mask, var='biomass', spacing=4000):
# fit_score = score(test[1][0], y_pred)

# make the grid
grid = chain.grid(spacing=spacing, region=region, data_names=[var], dims=('y', 'x'))
grid = chain.grid(spacing=spacing, region=region, data_names=[var], dims=("y", "x"))
grid = vd.distance_mask(
proj_coords,
maxdist=4 * spacing,
Expand Down
81 changes: 44 additions & 37 deletions carbonplan_forest_risks/load/cmip.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,74 +12,79 @@
from .. import setup, utils

members = {
'CanESM5-CanOE': 'r3i1p2f1',
'MIROC-ES2L': 'r1i1p1f2',
'ACCESS-CM2': 'r1i1p1f1',
'ACCESS-ESM1-5': 'r10i1p1f1',
'MRI-ESM2-0': 'r1i1p1f1',
'MPI-ESM1-2-LR': 'r10i1p1f1',
"CanESM5-CanOE": "r3i1p2f1",
"MIROC-ES2L": "r1i1p1f2",
"ACCESS-CM2": "r1i1p1f1",
"ACCESS-ESM1-5": "r10i1p1f1",
"MRI-ESM2-0": "r1i1p1f1",
"MPI-ESM1-2-LR": "r10i1p1f1",
}


@retry(stop=stop_after_attempt(5))
def cmip(
store='az',
store="az",
df=None,
tlim=None,
model=None,
scenario=None,
coarsen=None,
variables=['ppt', 'tmean'],
variables=["ppt", "tmean"],
mask=None,
member=None,
method='bias-corrected',
sampling='annual',
method="bias-corrected",
sampling="annual",
historical=False,
remove_nans=False,
):

with warnings.catch_warnings():
warnings.simplefilter('ignore', category=ResourceWarning)
warnings.simplefilter('ignore', category=FutureWarning)
warnings.simplefilter('ignore', category=RuntimeWarning)
warnings.simplefilter("ignore", category=ResourceWarning)
warnings.simplefilter("ignore", category=FutureWarning)
warnings.simplefilter("ignore", category=RuntimeWarning)

if scenario is None:
raise ValueError('must specify scenario')
raise ValueError("must specify scenario")
if model is None:
raise ValueError('must specify model')
raise ValueError("must specify model")
if member is None:
member = members[model]


path = setup.loading(store)

prefix = f'cmip6/{method}/conus/4000m/{sampling}/{model}.{scenario}.{member}.zarr'
prefix = (
f"cmip6/{method}/conus/4000m/{sampling}/{model}.{scenario}.{member}.zarr"
)

if store == 'az':
if store == "az":
mapper = zarr.storage.ABSStore(
'carbonplan-downscaling', prefix=prefix, account_name='carbonplan'
"carbonplan-downscaling", prefix=prefix, account_name="carbonplan"
)
else:
mapper = fsspec.get_mapper((path / 'carbonplan-downscaling' / prefix).as_uri())
mapper = fsspec.get_mapper(
(path / "carbonplan-downscaling" / prefix).as_uri()
)

ds = xr.open_zarr(mapper, consolidated=True)

if historical:
prefix = f'cmip6/{method}/conus/4000m/{sampling}/{model}.historical.{member}.zarr'
prefix = f"cmip6/{method}/conus/4000m/{sampling}/{model}.historical.{member}.zarr"

if store == 'az':
if store == "az":
mapper = zarr.storage.ABSStore(
'carbonplan-downscaling', prefix=prefix, account_name='carbonplan'
"carbonplan-downscaling", prefix=prefix, account_name="carbonplan"
)
else:
mapper = fsspec.get_mapper((path / 'carbonplan-downscaling' / prefix).as_uri())
mapper = fsspec.get_mapper(
(path / "carbonplan-downscaling" / prefix).as_uri()
)

ds_historical = xr.open_zarr(mapper, consolidated=True)

ds = xr.concat([ds_historical, ds], 'time')
ds = xr.concat([ds_historical, ds], "time")

ds['cwd'] = ds['def']
ds['pdsi'] = ds['pdsi'].clip(-16, 16)
ds["cwd"] = ds["def"]
ds["pdsi"] = ds["pdsi"].clip(-16, 16)

X = xr.Dataset()
keys = variables
Expand All @@ -99,29 +104,31 @@ def cmip(
if coarsen:
X_coarse = xr.Dataset()
for key in keys:
X_coarse[key] = X[key].coarsen(x=coarsen, y=coarsen, boundary='trim').mean()
X_coarse[key] = (
X[key].coarsen(x=coarsen, y=coarsen, boundary="trim").mean()
)
X = X_coarse

if df is not None:
t = Affine(*utils.albers_conus_transform(4000))
p1 = Proj(utils.albers_conus_crs())
p2 = Proj(proj='latlong', datum='WGS84')
x, y = transform(p2, p1, df['lon'].values, df['lat'].values)
p2 = Proj(proj="latlong", datum="WGS84")
x, y = transform(p2, p1, df["lon"].values, df["lat"].values)
rc = rowcol(t, x, y)
ind_r = xr.DataArray(rc[0], dims=['c'])
ind_c = xr.DataArray(rc[1], dims=['c'])
ind_r = xr.DataArray(rc[0], dims=["c"])
ind_c = xr.DataArray(rc[1], dims=["c"])

base = X[keys].isel(y=ind_r, x=ind_c).load()
for key in keys:
df[key + '_mean'] = base[key].mean('time').values
df[key + '_min'] = base[key].min('time').values
df[key + '_max'] = base[key].max('time').values
df[key + "_mean"] = base[key].mean("time").values
df[key + "_min"] = base[key].min("time").values
df[key + "_max"] = base[key].max("time").values
if remove_nans:
for key in keys:
df = df[~np.isnan(df[key + '_mean'])]
df = df[~np.isnan(df[key + "_mean"])]
df = df.reset_index(drop=True)
return df

X = X.drop(['x', 'y'])
X = X.drop(["x", "y"])
X.load(retries=10)
return X
Loading

0 comments on commit 1cd21d3

Please sign in to comment.