diff --git a/fpcup/_typing.py b/fpcup/_typing.py index a71cfdb..da0a72b 100644 --- a/fpcup/_typing.py +++ b/fpcup/_typing.py @@ -5,6 +5,12 @@ from os import PathLike from typing import Callable, Iterable, Optional, Type +from geopandas import GeoSeries +from shapely import Point, Polygon + Coordinates = tuple[RealNumber, RealNumber] PathOrStr = PathLike | str StringDict = dict[str, str] + +AreaDict = dict[str, Polygon] +BoundaryDict = dict[str, GeoSeries] diff --git a/fpcup/plotting.py b/fpcup/plotting.py index 274a2af..7b4a30d 100644 --- a/fpcup/plotting.py +++ b/fpcup/plotting.py @@ -10,47 +10,63 @@ from matplotlib import pyplot as plt, dates as mdates, patches as mpatches, ticker as mticker from matplotlib import colormaps, rcParams -rcParams.update({"axes.grid": True, "figure.dpi": 600, "grid.linestyle": "--", "hist.bins": 15, "legend.edgecolor": "black", "legend.framealpha": 1, "savefig.dpi": 600}) +rcParams.update({"axes.grid": True, + "figure.dpi": 600, + "grid.linestyle": "--", + "hist.bins": 15, + "image.cmap": "cividis", + "legend.edgecolor": "black", "legend.framealpha": 1, + "savefig.dpi": 600}) from mpl_toolkits.axes_grid1 import make_axes_locatable cividis_discrete = colormaps["cividis"].resampled(10) from ._brp_dictionary import brp_categories_colours, brp_crops_colours from ._typing import Iterable, Optional, PathOrStr, RealNumber, StringDict +from .analysis import KEYS_AGGREGATE +from .constants import CRS_AMERSFOORT, WGS84 from .model import Summary, parameter_names -from .province import nl_boundary, province_area, province_boundary, province_coarse - -# @mticker.FuncFormatter -# def _capitalise_ticks(x, pos): -# """ -# Helper function to capitalise str ticks in plots. -# """ -# try: -# return x.capitalize() -# except AttributeError: -# return "abc" -# _capitalise_ticks = mticker.StrMethodFormatter("abc{x:}") - -def plot_outline(ax: plt.Axes, province: str="All", *, coarse=False, **kwargs) -> None: +from .province import PROVINCE_NAMES, area, area_coarse, boundary, boundary_coarse, aggregate_h3 +from .tools import make_iterable + + +def plot_outline(ax: plt.Axes, province: str="Netherlands", *, + coarse: bool=False, crs: str=CRS_AMERSFOORT, **kwargs) -> None: """ - Plot an outline of the Netherlands ("All") or a specific province (e.g. "Zuid-Holland"). + Plot an outline of the Netherlands or a specific province (e.g. "Zuid-Holland"). """ - if province == "All": - boundary = nl_boundary + if coarse: + line = boundary_coarse[province].to_crs(crs) else: - if coarse: - boundary = gpd.GeoSeries(province_coarse[province].boundary) - else: - boundary = province_boundary[province] + line = boundary[province].to_crs(crs) line_kw = {"color": "black", "lw": 1, **kwargs} - boundary.plot(ax=ax, **line_kw) + line.plot(ax=ax, **line_kw) + + +def _configure_map_panels(axs: plt.Axes | Iterable[plt.Axes], province: str | Iterable[str]="Netherlands", **kwargs) -> None: + """ + Apply default settings to map panels. + **kwargs are passed to `plot_outline` - e.g. coarse, crs, ... + """ + axs = make_iterable(axs) + provinces = make_iterable(province) + for ax in axs: + # Country/Province outline(s) + for p in provinces: + plot_outline(ax, p, **kwargs) + + # Axis settings + ax.set_axis_off() + ax.axis("equal") + -def column_to_title(column: str) -> str: +def _column_to_title(column: str) -> str: """ Clean up a column name (e.g. "crop_species") so it can be used as a title (e.g. "Crop species"). """ return column.capitalize().replace("_", " ") + def brp_histogram(data: gpd.GeoDataFrame, column: str, *, figsize=(3, 5), usexticks=True, xlabel: Optional[str]="Crop", title: Optional[str]=None, top5=True, saveto: Optional[PathOrStr]=None, **kwargs) -> None: """ @@ -97,16 +113,16 @@ def brp_histogram(data: gpd.GeoDataFrame, column: str, *, plt.show() plt.close() + def brp_map(data: gpd.GeoDataFrame, column: str, *, - province: Optional[str]="All", figsize=(10, 10), title: Optional[str]=None, rasterized=True, colour_dict: Optional[StringDict]=None, saveto: Optional[PathOrStr]=None, **kwargs) -> None: + province: str="Netherlands", figsize=(10, 10), title: Optional[str]=None, rasterized=True, colour_dict: Optional[StringDict]=None, saveto: Optional[PathOrStr]=None, **kwargs) -> None: """ Create a map of BRP polygons in the given column. If `province` is provided, only data within that province will be plotted, with the corresponding outline. """ # Select province data if desired - if province != "All": + if province != "Netherlands": assert "province" in data.columns, f"Cannot plot data by province - data do not have a 'province' column\n(columns: {data.columns}" - province = province.title() data = data.loc[data["province"] == province] # Create figure @@ -119,24 +135,22 @@ def brp_map(data: gpd.GeoDataFrame, column: str, *, # Generate dummy patches with the same colour mapping and add those to the legend colour_patches = [mpatches.Patch(color=colour, label=label.capitalize()) for label, colour in colour_dict.items() if label in data[column].unique()] - ax.legend(handles=colour_patches, loc="lower right", fontsize=12, title=column_to_title(column)) + ax.legend(handles=colour_patches, loc="lower right", fontsize=12, title=_column_to_title(column)) # If colours are not specified, simply plot the data and let geopandas handle the colours else: data.plot(ax=ax, column=column, rasterized=rasterized, **kwargs) - # Add a country/province outline - plot_outline(ax, province) - + # Panel settings + _configure_map_panels(ax, province) ax.set_title(title) - ax.set_axis_off() - ax.axis("equal") if saveto: plt.savefig(saveto, bbox_inches="tight") plt.show() plt.close() + def brp_crop_map_split(data: gpd.GeoDataFrame, column: str="crop_species", *, crops: Iterable[str]=brp_crops_colours.keys(), figsize=(14, 3.5), shape=(1, 5), title: Optional[str]=None, rasterized=True, saveto: Optional[PathOrStr]=None, **kwargs) -> None: """ @@ -151,12 +165,9 @@ def brp_crop_map_split(data: gpd.GeoDataFrame, column: str="crop_species", *, data_here = data.loc[data[column] == crop] data_here.plot(ax=ax, color="black", rasterized=rasterized, **kwargs) - if nl_boundary is not None: - nl_boundary.plot(ax=ax, color="black", lw=0.2) - ax.set_title(crop.capitalize(), color=colour) - ax.set_axis_off() - ax.axis("equal") + + _configure_map_panels(axs.ravel(), lw=0.2) fig.suptitle(title) @@ -165,6 +176,7 @@ def brp_crop_map_split(data: gpd.GeoDataFrame, column: str="crop_species", *, plt.show() plt.close() + def replace_year_in_datetime(date: dt.date, newyear: int=2000) -> dt.date: """ For a datetime object yyyy-mm-dd, replace yyyy with newyear. @@ -172,6 +184,7 @@ def replace_year_in_datetime(date: dt.date, newyear: int=2000) -> dt.date: """ return date.replace(year=newyear) + def plot_wofost_ensemble_results(outputs: Iterable[pd.DataFrame], keys: Iterable[str]=None, *, title: Optional[str]=None, saveto: Optional[PathOrStr]=None, replace_years=True, progressbar=True, leave_progressbar=False) -> None: """ @@ -214,6 +227,7 @@ def plot_wofost_ensemble_results(outputs: Iterable[pd.DataFrame], keys: Iterable plt.close() + def _numerical_or_date_bins(column: pd.Series) -> int | pd.DatetimeIndex: """ Generate bins for a column based on its data type. @@ -223,15 +237,12 @@ def _numerical_or_date_bins(column: pd.Series) -> int | pd.DatetimeIndex: else: return rcParams["hist.bins"] -def plot_wofost_ensemble_summary(summary: Summary, keys: Iterable[str]=None, *, - weights: Optional[Iterable[RealNumber]]=None, title: Optional[str]=None, province: Optional[str]="All", saveto: Optional[PathOrStr]=None) -> None: + +def plot_wofost_ensemble_summary(summary: Summary, keys: Iterable[str]=KEYS_AGGREGATE, *, + aggregate: bool=True, weights: Optional[Iterable[RealNumber]]=None, title: Optional[str]=None, province: Optional[str]="Netherlands", saveto: Optional[PathOrStr]=None) -> None: """ - Plot WOFOST ensemble results. + Plot histograms and (aggregate) maps showing WOFOST run summaries. """ - # If no keys were specified, get all of them - if keys is None: - keys = summary.keys() - # Create figure and panels fig, axs = plt.subplots(nrows=2, ncols=len(keys), sharey="row", figsize=(15, 5), gridspec_kw={"hspace": 0.25, "height_ratios": [1, 1.5]}) @@ -267,12 +278,7 @@ def plot_wofost_ensemble_summary(summary: Summary, keys: Iterable[str]=None, *, im = summary.plot(column, ax=ax_col[1], rasterized=True, vmin=vmin_here, vmax=vmax_here, legend=True, cax=cax, cmap=cividis_discrete, legend_kwds={"location": "bottom"}) # Settings for map panels - for ax in axs[1]: - # Add a country/province outline - plot_outline(ax, province) - - ax.set_axis_off() - ax.axis("equal") + _configure_map_panels(axs[1]) axs[0, 0].set_ylabel("Distribution") fig.align_xlabels() @@ -321,13 +327,7 @@ def plot_wofost_ensemble_summary_aggregate(aggregate: gpd.GeoDataFrame, keys: It im = aggregate.plot(column, ax=ax, rasterized=True, vmin=vmin_here, vmax=vmax_here, legend=True, cax=cax, cmap=cividis_discrete, legend_kwds={"location": "bottom", "label": key}) # Settings for map panels - for ax in axs: - # Add a country/province outline - for province in province_coarse.keys(): - plot_outline(ax, province, coarse=True, lw=0.5) - - ax.set_axis_off() - ax.axis("equal") + _configure_map_panels(axs, PROVINCE_NAMES, coarse=True, lw=0.5) if title is None: title = f"Results aggregated by province" diff --git a/fpcup/province.py b/fpcup/province.py index 71efd45..c09c28a 100644 --- a/fpcup/province.py +++ b/fpcup/province.py @@ -1,33 +1,46 @@ """ -(Try to) load map backgrounds from file so they can be plotted. +Geography-related constants and methods. +Includes polygons/outlines of the Netherlands and its provinces. """ - import geopandas as gpd +gpd.options.io_engine = "pyogrio" +import h3pandas import numpy as np from pandas import DataFrame, Series from tqdm import tqdm -from ._typing import Iterable, Optional -from .constants import CRS_AMERSFOORT +from ._typing import AreaDict, BoundaryDict, Callable, Iterable, Optional +from .constants import CRS_AMERSFOORT, WGS84 from .settings import DEFAULT_DATA -# Load the outline of the Netherlands -nl = gpd.read_file(DEFAULT_DATA/"NL_borders.geojson") -nl_boundary = nl.boundary +# Constants +PROVINCE_NAMES = ("Fryslân", "Gelderland", "Noord-Brabant", "Noord-Holland", "Overijssel", "Zuid-Holland", "Groningen", "Zeeland", "Drenthe", "Flevoland", "Limburg", "Utrecht") # Sorted by area +NETHERLANDS = "Netherlands" +NAMES = PROVINCE_NAMES + (NETHERLANDS, ) + +# Load the Netherlands shapefile +_netherlands = gpd.read_file(DEFAULT_DATA/"NL_borders.geojson") +_area_netherlands = {NETHERLANDS: _netherlands.iloc[0].geometry} # Polygon - to be used in comparisons +_boundary_netherlands = {NETHERLANDS: _netherlands.boundary} # GeoSeries - to be used with .plot() -# Load the provinces +# Load the province shapefiles _provinces = gpd.read_file(DEFAULT_DATA/"NL_provinces.geojson") _provinces_coarse = gpd.read_file(DEFAULT_DATA/"NL_provinces_coarse.geojson") -province_names = list(_provinces["naamOfficieel"]) + ["Friesland"] -# Access individual provinces using a dictionary, e.g. province_boundary["Zuid-Holland"] -province_area = {name: poly for name, poly in zip(_provinces["naamOfficieel"], _provinces["geometry"])} -province_boundary = {name: gpd.GeoSeries(outline) for name, outline in zip(_provinces["naamOfficieel"], _provinces.boundary)} -province_coarse = {name: poly for name, poly in zip(_provinces_coarse["naamOfficieel"], _provinces_coarse["geometry"])} +# Access individual provinces using a dictionary, e.g. area["Zuid-Holland"] +# Note: these contain bare Polygon/MultiPolygon objects, with no CRS. +area = {name: poly for name, poly in zip(_provinces["naamOfficieel"], _provinces["geometry"])} +area = {**_area_netherlands, **area} + +area_coarse = {name: poly for name, poly in zip(_provinces_coarse["naamOfficieel"], _provinces_coarse["geometry"])} + +# Access individual provinces using a dictionary, e.g. boundary["Zuid-Holland"] +# Note: these contain GeoSeries objects with 1 entry, with a CRS set. +boundary = {name: gpd.GeoSeries(outline, crs=CRS_AMERSFOORT) for name, outline in zip(_provinces["naamOfficieel"], _provinces.boundary)} +boundary = {**_boundary_netherlands, **boundary} + +boundary_coarse = {name: gpd.GeoSeries(poly.boundary, crs=CRS_AMERSFOORT) for name, poly in area_coarse.items()} -# Add an alias for Friesland/Fryslân -province_area["Friesland"] = province_area["Fryslân"] -province_boundary["Friesland"] = province_boundary["Fryslân"] def process_input_province(province: str) -> str: """ @@ -39,32 +52,40 @@ def process_input_province(province: str) -> str: # Apply common alternatives if province in ("Friesland", "Fryslan"): province = "Fryslân" + elif province in ("the Netherlands", "NL", "All"): + province = "Netherlands" return province -def is_in_province(data: gpd.GeoDataFrame, province: str, *, province_data: dict=province_coarse, use_centroid=True) -> Iterable[bool]: + +def is_in_province(_data: gpd.GeoDataFrame, province: str, *, + province_data: AreaDict=area_coarse, use_centroid=True) -> Iterable[bool]: """ For a series of geometries (e.g. BRP plots), determine if they are in the given province. Enable `use_centroid` to see if the centre of each plot falls within the province rather than the entire plot - this is useful for plots that are split between provinces. """ + assert province in PROVINCE_NAMES, f"Unknown province '{province}'." + area = province_data[province] if use_centroid: - data_here = data.centroid + data = _data.centroid else: - data_here = data + data = _data # Step 1: use the convex hull for a coarse selection - selection_coarse = data_here.within(area.convex_hull) + selection_coarse = data.within(area.convex_hull) # Step 2: use the real shape of the province - selection_fine = data_here.loc[selection_coarse].within(area) + selection_fine = data.loc[selection_coarse].within(area) # Update the coarse selection with the new information selection_coarse.loc[selection_coarse] = selection_fine return selection_coarse -def add_provinces(data: gpd.GeoDataFrame, *, new_column: str="province", province_data: dict=province_coarse, progressbar=True, leave_progressbar=True, **kwargs) -> None: + +def add_provinces(data: gpd.GeoDataFrame, *, + new_column: str="province", province_data: AreaDict=area_coarse, progressbar=True, leave_progressbar=True, **kwargs) -> None: """ Add a column with province names. Note: can get very slow for long dataframes. @@ -73,7 +94,8 @@ def add_provinces(data: gpd.GeoDataFrame, *, new_column: str="province", provinc data_crs = data.to_crs(CRS_AMERSFOORT) # Generate an empty Series which will be populated with time - province_list = Series(data=np.tile("", len(data_crs)), name="province", dtype=str, index=data_crs.index) + province_list = Series(data=np.tile("", len(data_crs)), + name="province", dtype=str, index=data_crs.index) # Loop over the provinces, find the relevant entries, and fill in the list for province_name in tqdm(province_data.keys(), desc="Assigning labels", unit="province", disable=not progressbar, leave=leave_progressbar): @@ -93,7 +115,9 @@ def add_provinces(data: gpd.GeoDataFrame, *, new_column: str="province", provinc # Add the series to the dataframe data[new_column] = province_list -def add_province_geometry(data: DataFrame, which: str="area", *, column_name: Optional[str]=None, crs: str=CRS_AMERSFOORT) -> gpd.GeoDataFrame: + +def add_province_geometry(data: DataFrame, which: str="area", *, + column_name: Optional[str]=None, crs: str=CRS_AMERSFOORT) -> gpd.GeoDataFrame: """ Add a column with province geometry to a DataFrame with a province name column/index. """ @@ -107,7 +131,7 @@ def add_province_geometry(data: DataFrame, which: str="area", *, column_name: Op column = data[column_name] # Apply the dictionary mapping - geometry = column.map(province_area) + geometry = column.map(area) data_new = gpd.GeoDataFrame(data, geometry=geometry, crs=crs) # Change to outline if desired @@ -120,3 +144,32 @@ def add_province_geometry(data: DataFrame, which: str="area", *, column_name: Op raise ValueError(f"Cannot add geometries of type `{which}`") return data_new + + +def aggregate_h3(_data: gpd.GeoDataFrame, functions: dict[str, Callable] | Callable | str="mean", *, + level: int=6, clipto: Optional[str]="Netherlands") -> gpd.GeoDataFrame: + """ + Aggregate data to the H3 hexagonal grid. + `functions` is passed to DataFrame.agg. + `clipto` is used to get a geometry, e.g. the Netherlands or one province, to clip the results to. Set it to `None` to preserve the full grid. + + TO DO: See if using `clipto` to filter data to a given province before aggregation is faster. + """ + # Convert the input to WGS84 for the aggregation + crs_original = _data.crs + data = _data.copy() + data["geometry"] = data.centroid.to_crs(WGS84) + + # Aggregate the data and convert back to the original CRS + data_h3 = data.h3.geo_to_h3_aggregate(level, functions) + data_h3.to_crs(crs_original, inplace=True) + + # Clip the data if desired + if clipto is not None: + assert clipto in area.keys(), f"Cannot clip the H3 grid to '{clipto}' - unknown name. Please provide a province name, 'Netherlands', or None." + clipto_geometry = area[clipto] + + data_h3["geometry"] = data_h3.intersection(clipto_geometry) + data_h3 = data_h3.loc[~data_h3.is_empty] + + return data_h3 diff --git a/plot_results.py b/plot_results.py index 2647d9d..1b99376 100644 --- a/plot_results.py +++ b/plot_results.py @@ -11,12 +11,12 @@ parser.add_argument("output_dir", help="folder to load PCSE outputs from", type=Path) parser.add_argument("-y", "--replace_years", help="replace all years in the output with 2000", action="store_true") parser.add_argument("--vector_max", help="number of runs at which to switch from vector (PDF) to raster (PNG) files", type=int, default=5000) -parser.add_argument("-p", "--province", help="province to select plots from (or all)", default="All", choices=fpcup.province.province_names+["All"], type=fpcup.province.process_input_province) +parser.add_argument("-p", "--province", help="province to select plots from (or all of the Netherlands)", default="Netherlands", choices=fpcup.province.NAMES, type=fpcup.province.process_input_province) parser.add_argument("-s", "--sample", help="load only a subsample of the outputs, for testing", action="store_true") parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") args = parser.parse_args() -SINGLE_PROVINCE = (args.province != "All") +SINGLE_PROVINCE = (args.province != "Netherlands") # Set up the input/output directories tag = args.output_dir.stem @@ -36,12 +36,6 @@ if args.verbose: print(f"Loaded summary file -- {len(summary)} rows") -# Check if area information is available (to be used in weights later) -AREA_AVAILABLE = ("area" in summary.columns) -if args.verbose: - un = "un" if not AREA_AVAILABLE else "" - print(f"Plot areas {un}available; histograms and means will be {un}weighted") - # Add province information if this is not available if "province" not in summary.columns: fpcup.province.add_provinces(summary, leave_progressbar=args.verbose) @@ -55,18 +49,27 @@ if args.verbose: print(f"Selected only sites in {args.province} -- {len(summary)} rows") -# Use areas for weights +# Check if area information is available (to be used in weights) +AREA_AVAILABLE = ("area" in summary.columns) areas = summary["area"] if AREA_AVAILABLE else None +if args.verbose: + un = "un" if not AREA_AVAILABLE else "" + print(f"Plot areas {un}available; histograms and means will be {un}weighted") + +# Define aggregation functions +if AREA_AVAILABLE: + aggregator = fpcup.analysis.weighted_mean_dict(summary) +else: + aggregator = {key: "mean" for key in fpcup.analysis.KEYS_AGGREGATE} # Useful information if args.verbose: print(f"Figures will be saved as -{tag}") # Plot summary results -keys_to_plot = ["LAIMAX", "TWSO", "CTRAT", "CEVST", "DOE", "DOM"] filename_summary = results_dir / f"WOFOST_{tag}-summary.pdf" -fpcup.plotting.plot_wofost_ensemble_summary(summary, keys=keys_to_plot, weights=areas, saveto=filename_summary, title=f"Summary of {len(summary)} WOFOST runs: {tag}", province=args.province) +fpcup.plotting.plot_wofost_ensemble_summary(summary, weights=areas, saveto=filename_summary, title=f"Summary of {len(summary)} WOFOST runs: {tag}", province=args.province) if args.verbose: print(f"Saved batch results plot to {filename_summary.absolute()}") @@ -75,19 +78,17 @@ # Calculate the mean per province of several variables, weighted by plot area if possible if AREA_AVAILABLE: - mean_func = fpcup.analysis.weighted_mean_dict(summary) filename_means = results_dir / f"WOFOST_{tag}-weighted-mean.csv" if args.verbose: print("Calculating weighted means") # Use a normal mean if there is no area information available else: - mean_func = {key: "mean" for key in fpcup.analysis.KEYS_AGGREGATE} filename_means = results_dir / f"WOFOST_{tag}-mean.csv" if args.verbose: print("Could not calculate weighted means because there is no 'area' column -- defaulting to a regular mean") -byprovince_mean = byprovince[fpcup.analysis.KEYS_AGGREGATE].agg(mean_func) +byprovince_mean = byprovince[fpcup.analysis.KEYS_AGGREGATE].agg(aggregator) byprovince_mean.to_csv(filename_means) if args.verbose: print(f"Saved aggregate mean columns to {filename_means.absolute()}") @@ -100,20 +101,6 @@ if args.verbose: print(f"Saved aggregate mean plot to {filename_aggregate.absolute()}") -# Use H3 for aggregation -import h3pandas -summary2 = summary.copy() -summary2["geometry"] = summary.centroid.to_crs(fpcup.constants.WGS84) -summaryh3 = summary2.h3.geo_to_h3_aggregate(5, mean_func) -summaryh3.to_crs(fpcup.constants.CRS_AMERSFOORT, inplace=True) - -from matplotlib import pyplot as plt -fig, ax = plt.subplots(figsize=(8,8)) -summaryh3.plot("LAIMAX", ax=ax) -fpcup.plotting.plot_outline(ax) -plt.savefig("test.pdf") -plt.close() - # Space between summary and outputs sections if args.verbose: print() diff --git a/process_brp.py b/process_brp.py index 9e4bae2..8575cb8 100644 --- a/process_brp.py +++ b/process_brp.py @@ -6,6 +6,7 @@ %run process_brp.py data/brp/brpgewaspercelen_definitief_2022.gpkg -vp """ import geopandas as gpd +gpd.options.io_engine = "pyogrio" gpd.pd.options.mode.chained_assignment = None # Prevents unneeded warnings import fpcup diff --git a/process_mars.py b/process_mars.py index d85de72..b276f9a 100644 --- a/process_mars.py +++ b/process_mars.py @@ -19,7 +19,7 @@ from pcse.util import wind10to2 -from fpcup.province import nl, CRS_AMERSFOORT +import fpcup filename = Path("data/meteo/NL-weather_ver2024_01_28946_162168065.csv") @@ -42,7 +42,7 @@ data = gpd.GeoDataFrame(data, geometry=gpd.points_from_xy(data["LONGITUDE"], data["LATITUDE"], crs="WGS84")) # Regrid: definitions -data.to_crs(CRS_AMERSFOORT, inplace=True) +data.to_crs(fpcup.constants.CRS_AMERSFOORT, inplace=True) xmin, xmax, ymin, ymax = data.total_bounds cellsize = 250 # m newx = np.arange(xmin, xmax+cellsize, cellsize) @@ -63,7 +63,7 @@ fig, ax = plt.subplots(1, 1, figsize=(5,5)) im = ax.imshow(gridded, origin="lower", extent=(newx.min(), newx.max(), newy.min(), newy.max()), vmin=vmin, vmax=vmax) data_day.plot(key, ax=ax, edgecolor="black", vmin=vmin, vmax=vmax) -nl.boundary.plot(ax=ax, color="black", lw=1) +fpcup.plotting.plot_outline(ax) fig.colorbar(im, ax=ax, label=key) diff --git a/requirements.txt b/requirements.txt index 553e0c2..b128bea 100644 --- a/requirements.txt +++ b/requirements.txt @@ -3,3 +3,4 @@ pandas>=2.1.3 pcse>=5.5.5 tqdm>=4.66.1 geopandas>=0.14.2 +h3pandas>=0.2.6 \ No newline at end of file diff --git a/wofost_brp.py b/wofost_brp.py index b0dbab6..0853674 100644 --- a/wofost_brp.py +++ b/wofost_brp.py @@ -11,13 +11,15 @@ parser = argparse.ArgumentParser(description="Run PCSE for plots within the BRP.") parser.add_argument("brp_filename", help="file to load the BRP from", type=fpcup.io.Path) parser.add_argument("-c", "--crop", help="crop to run simulations on (or all)", default="All", choices=("barley", "maize", "sorghum", "soy", "wheat"), type=str.lower) -parser.add_argument("-p", "--province", help="province to select plots from (or all)", default="All", choices=fpcup.province.province_names+["All"], type=fpcup.province.process_input_province) +parser.add_argument("-p", "--province", help="province to select plots from (or all)", default="Netherlands", choices=fpcup.province.NAMES, type=fpcup.province.process_input_province) parser.add_argument("-d", "--data_dir", help="folder to load PCSE data from", type=fpcup.io.Path, default=fpcup.settings.DEFAULT_DATA) parser.add_argument("-o", "--output_dir", help="folder to save PCSE outputs to", type=fpcup.io.Path, default=None) parser.add_argument("-f", "--force", help="run all models even if the associated file already exists", action="store_true") parser.add_argument("-v", "--verbose", help="increase output verbosity", action="store_true") args = parser.parse_args() +SINGLE_PROVINCE = (args.province != "Netherlands") + # Get the year from the BRP filename year = int(args.brp_filename.stem.split("_")[-1].split("-")[0]) @@ -40,7 +42,7 @@ print(f"Loaded BRP data from {args.brp_filename.absolute()} -- {len(brp)} sites") # If we are only doing one province, select only the relevant lines from the BRP file -if args.province != "All": +if SINGLE_PROVINCE: brp = brp.loc[brp["province"] == args.province] if args.verbose: