diff --git a/src/natcap/invest/carbon/reporter.py b/src/natcap/invest/carbon/reporter.py index 115c773e09..9321786a88 100644 --- a/src/natcap/invest/carbon/reporter.py +++ b/src/natcap/invest/carbon/reporter.py @@ -208,10 +208,6 @@ def report(file_registry: dict, args_dict: dict, model_spec: ModelSpec, agg_results_table = _generate_agg_results_table(args_dict, file_registry) - lulc_pre_caption = gettext( - 'Values in the legend are listed in order of frequency (most common ' - 'first).') - with open(target_html_filepath, 'w', encoding='utf-8') as target_file: target_file.write(TEMPLATE.render( report_script=model_spec.reporter, @@ -226,7 +222,7 @@ def report(file_registry: dict, args_dict: dict, model_spec: ModelSpec, agg_results_table=agg_results_table, inputs_img_src=inputs_img_src, inputs_caption=input_raster_caption, - lulc_pre_caption=lulc_pre_caption, + lulc_pre_caption=report_constants.LULC_PRE_CAPTION, outputs_img_src=outputs_img_src, outputs_caption=output_raster_caption, intermediate_raster_sections=intermediate_raster_sections, diff --git a/src/natcap/invest/reports/raster_utils.py b/src/natcap/invest/reports/raster_utils.py index 3bb8152d92..202caea4e8 100644 --- a/src/natcap/invest/reports/raster_utils.py +++ b/src/natcap/invest/reports/raster_utils.py @@ -170,6 +170,10 @@ def __post_init__(self): self.colormap = plt.get_cmap(self.colormap if self.colormap else COLORMAPS[self.datatype]) + def plot_on_axis(self, fig, ax, arr, cmap, imshow_kwargs, colorbar_kwargs): + mappable = ax.imshow(arr, cmap=cmap, **imshow_kwargs) + fig.colorbar(mappable, ax=ax, **colorbar_kwargs) + def build_raster_plot_configs(id_lookup_table, raster_plot_tuples): """Build RasterPlotConfigs for use in plotting input or output rasters. @@ -462,8 +466,9 @@ def plot_raster_list(raster_list: list[RasterPlotConfig]): len(patches), n_plots, xy_ratio)) leg.set_in_layout(True) else: - mappable = ax.imshow(arr, cmap=cmap, **imshow_kwargs) - fig.colorbar(mappable, ax=ax, **colorbar_kwargs) + config.plot_on_axis( + fig, ax, arr, cmap, imshow_kwargs, colorbar_kwargs) + [ax.set_axis_off() for ax in axs.flatten()] return fig diff --git a/src/natcap/invest/reports/report_constants.py b/src/natcap/invest/reports/report_constants.py index d8fa3e31dd..1c4db454fd 100644 --- a/src/natcap/invest/reports/report_constants.py +++ b/src/natcap/invest/reports/report_constants.py @@ -18,4 +18,9 @@ 'in GIS to assess its accuracy.' ) +LULC_PRE_CAPTION = gettext( + 'LULC values in the legend are listed in order of frequency (most common ' + 'first).' +) + TABLE_PAGINATION_THRESHOLD = 10 diff --git a/src/natcap/invest/reports/templates/models/urban_mental_health.html b/src/natcap/invest/reports/templates/models/urban_mental_health.html new file mode 100644 index 0000000000..ca515772f1 --- /dev/null +++ b/src/natcap/invest/reports/templates/models/urban_mental_health.html @@ -0,0 +1,149 @@ +{% extends 'base.html' %} + +{% block styles %} + {{ super() }} + {% include 'datatable-styles.html' %} +{% endblock styles %} + +{% block content %} + + {{ super() }} + + {% from 'args-table.html' import args_table %} + {% from 'caption.html' import caption %} + {% from 'content-grid.html' import content_grid %} + {% from 'metadata.html' import list_metadata %} + {% from 'raster-plot-img.html' import raster_plot_img %} + {% from 'wide-table.html' import wide_table %} + + + +

Results

+ + {{ accordion_section( + 'Aggregate Results', + agg_results_table | safe + )}} + + {{ accordion_section( + 'Primary Outputs', + content_grid([ + (caption(raster_group_caption, pre_caption=True), 100), + (raster_plot_img(outputs_img_src, 'Primary Outputs'), 100), + (caption(outputs_caption, definition_list=True), 100) + ]) + )}} + + {{ accordion_section( + 'Preventable Cases by AOI Vector', + content_grid( + [ + (content_grid([ + ('
', 100), + (caption(cases_map_caption, aggregate_map_source_list), 100) + ]), 50), + + (content_grid([ + ('
', 100), + (caption(cost_map_caption, aggregate_map_source_list), 100) + ]), 50) + ] if cost_map_json else [ + (content_grid([ + ('
', 100), + (caption(cases_map_caption, aggregate_map_source_list), 100) + ]), 100) + ] + ) + ) }} + + {% for section in intermediate_raster_sections %} + {{ accordion_section( + section['heading'], + content_grid([ + (caption(raster_group_caption, pre_caption=True), 100), + (raster_plot_img(section['img_src'], section['heading']), 100), + (caption(section['caption'], definition_list=True), 100) + ]) + )}} + {% endfor %} + + {{ accordion_section( + 'Output Raster Stats', + content_grid([ + (stats_table_note, 100), + (wide_table( + output_raster_stats_table | safe, + font_size_px=16 + ), 100) + ]) + )}} + +

Inputs

+ + {% if args_dict != None %} + {{ accordion_section( + 'Arguments', + args_table(args_dict) + )}} + {% endif %} + + {% set input_items = [ + (caption(raster_group_caption, pre_caption=True), 100) + ] %} + + {% if args_dict['lulc_base'] %} + {% set input_items = input_items + [ + (caption(lulc_pre_caption, pre_caption=True), 100) + ] %} + {% endif %} + + {% set input_items = input_items + [ + (raster_plot_img(inputs_img_src, 'Raster Inputs'), 100), + (caption(inputs_caption, definition_list=True), 100) + ] %} + + {{ accordion_section( + 'Input Maps', + content_grid(input_items) + )}} + + {{ accordion_section( + 'Input Raster Stats', + content_grid([ + (stats_table_note, 100), + (wide_table( + input_raster_stats_table | safe, + font_size_px=16 + ), 100) + ]) + )}} + +

Metadata

+ + {{ + accordion_section( + 'Output Filenames and Descriptions', + list_metadata(model_spec_outputs), + expanded=False + ) + }} + +{% endblock content %} + + + +{% from 'vegalite-plot.html' import embed_vega %} + +{% block scripts %} + {{ super() }} + + + + {% include 'vega-embed-js.html' %} + + {% set chart_spec_id_list = [(cases_map_json, 'cases_map')] %} + {% if cost_map_json %} + {% set chart_spec_id_list = chart_spec_id_list + [(cost_map_json, 'cost_map')] %} + {% endif %} + {{ embed_vega(chart_spec_id_list) }} +{% endblock scripts %} diff --git a/src/natcap/invest/urban_mental_health/reporter.py b/src/natcap/invest/urban_mental_health/reporter.py new file mode 100644 index 0000000000..87c1462bc5 --- /dev/null +++ b/src/natcap/invest/urban_mental_health/reporter.py @@ -0,0 +1,500 @@ +import logging +import time + +import altair +import geopandas +import matplotlib +import numpy +import pandas +import pygeoprocessing +from pydantic import ConfigDict +from pydantic.dataclasses import dataclass +from osgeo import gdal + +from natcap.invest import __version__ +from natcap.invest import gettext +from natcap.invest.reports import jinja_env, raster_utils, report_constants, \ + vector_utils +from natcap.invest.spec import ModelSpec, FileRegistry + +from natcap.invest.reports.raster_utils import RasterDatatype, \ + RasterPlotConfig, RasterTransform + +LOGGER = logging.getLogger(__name__) + +TEMPLATE = jinja_env.get_template('models/urban_mental_health.html') + +MAP_WIDTH = 450 # pixels + +NEAR_ZERO_RANGE = (-0.01, 0.01) + +@dataclass(config=ConfigDict(arbitrary_types_allowed=True)) +class SpecialRangeRasterPlotConfig(RasterPlotConfig): + """RasterPlotConfig to allow for special handling of near-zero pixels.""" + special_value_range: tuple[float, float] | None = None + special_value_color: str | None = None + special_value_label: str | None = None + + def plot_on_axis(self, fig, ax, arr, cmap, imshow_kwargs, colorbar_kwargs): + if self.special_value_range is None: + LOGGER.info( + "No special value range for %s, using default plotting", + self.raster_path + ) + super().plot_on_axis( + fig, ax, arr, cmap, imshow_kwargs, colorbar_kwargs) + return + + low, high = self.special_value_range + + special_mask = ( + ~numpy.isnan(arr) & + (arr >= low) & + (arr <= high) + ) + + main_arr = numpy.array(arr, copy=True) + main_arr[special_mask] = numpy.nan + + mappable = ax.imshow(main_arr, cmap=cmap, **imshow_kwargs) + cbar = fig.colorbar(mappable, ax=ax, **colorbar_kwargs) + + special_arr = numpy.full(arr.shape, numpy.nan) + special_arr[special_mask] = 1 + + ax.imshow( + special_arr, + cmap=matplotlib.colors.ListedColormap([self.special_value_color]), + interpolation='none', + vmin=1, + vmax=1 + ) + + patch = matplotlib.patches.Patch( + color=self.special_value_color, + label=self.special_value_label or f'{low} - {high}' + ) + + cbar_ax = cbar.ax + cbar_ax.legend( + handles=[patch], + loc='upper center', + bbox_to_anchor=(0.5, -0.05), + frameon=False, + ncol=1 + ) + + +def infer_continuous_or_divergent(raster_path: str) -> str: + """Infer if raster should have a 'continuous' or 'divergent' color ramp. + + Rules: + - If min value < ~0 --> 'divergent' + - Else --> 'continuous' + + Args: + raster_path (str): Path to raster. + + Returns: + str: 'continuous' or 'divergent' + """ + raster_info = pygeoprocessing.get_raster_info(raster_path) + nodata = raster_info['nodata'][0] + ds = gdal.OpenEx(raster_path, gdal.OF_RASTER) + band = ds.GetRasterBand(1) + + stats = band.GetStatistics(False, True) # (approx_ok, force) + min_val = stats[0] + LOGGER.info("Stats for %s: min=%s, nodata=%s", raster_path, min_val, nodata) + + if min_val is None: + arr = pygeoprocessing.raster_to_numpy_array(raster_path) + if nodata is not None: + valid_mask = ~numpy.isclose(arr, nodata, equal_nan=True) + valid_values = arr[valid_mask] + else: + valid_values = arr[~numpy.isnan(arr)] + + if valid_values.size == 0: + LOGGER.warning(f"No valid pixels found in {raster_path}, " + "defaulting to continuous") + return RasterDatatype.continuous + + min_val = round(numpy.nanmin(valid_values), 4) + + return RasterDatatype.divergent if min_val < NEAR_ZERO_RANGE[0] else RasterDatatype.continuous + + +def _get_conditional_raster_plot_tuples(model_spec: ModelSpec, + args_dict: dict, + file_registry: FileRegistry) -> tuple[ + list[tuple[str, ...]], + list[tuple[str, ...]], + list[list[tuple[str, ...]]]]: + + """ + inputs + - population raster + - baseline prevalence vector + + if NDVI: + inputs: + - ndvi_bas_path + - ndvi_alt_path + # - if lulc_attr_csv: show lulc_attr_csv and lulc_base + # - if lulc_alt: show lulc_alt as well + + if LULC: + inputs: + - lulc_base + - lulc_alt + - if not lulc_attr_csv: show lulc_attr_table that was built + - if ndvi_base: show it as well + + intermediates: + - LULC baseline reclassified to NDVI + - if lulc_alt: LULC alternate reclassified to NDVI + + + intermediates: + - ndvi_delta + - baseline cases + - baseline prevalence + + outputs: + - preventable_cases + - preventable_cost + - preventable cases cost sum table (just total cases and total cost) + - vector symbolized with field 'sum_cost' + - vector symbolized with field 'sum_cases' + + + """ + ndvi_colorramp = 'viridis_r' + + input_raster_config_list = [] + intermediate_raster_config_lists = [[ + RasterPlotConfig( + raster_path=file_registry['ndvi_base_buffer_mean_clipped'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_output('ndvi_base_buffer_mean_clipped'), + colormap=ndvi_colorramp + ), + RasterPlotConfig( + raster_path=file_registry['ndvi_alt_buffer_mean_clipped'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_output('ndvi_alt_buffer_mean_clipped'), + colormap=ndvi_colorramp + ), + RasterPlotConfig( + raster_path=file_registry['delta_ndvi'], + datatype=RasterDatatype.divergent, + spec=model_spec.get_output('delta_ndvi'), + )], + [RasterPlotConfig( + raster_path=file_registry['baseline_prevalence_raster'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_output('baseline_prevalence_raster') + ), + RasterPlotConfig( + raster_path=file_registry['baseline_cases'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_output('baseline_cases') + ) + ]] + + if args_dict['lulc_base']: + input_raster_config_list.append( + RasterPlotConfig( + raster_path=args_dict['lulc_base'], + datatype=RasterDatatype.nominal, + spec=model_spec.get_input('lulc_base') + ) + ) + if args_dict['lulc_alt']: + input_raster_config_list.append( + RasterPlotConfig( + raster_path=args_dict['lulc_alt'], + datatype=RasterDatatype.nominal, + spec=model_spec.get_input('lulc_alt') + ) + ) + + if args_dict["scenario"] == 'ndvi': + input_raster_config_list.insert(0, + RasterPlotConfig( + raster_path=args_dict['ndvi_base'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_input('ndvi_base'), + colormap=ndvi_colorramp + ) + ) + input_raster_config_list.insert(1, + RasterPlotConfig( + raster_path=args_dict['ndvi_alt'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_input('ndvi_alt'), + colormap=ndvi_colorramp + ) + ) + + else: + intermediate_raster_config_lists.insert(0, [ + RasterPlotConfig( + raster_path=file_registry['ndvi_base_aligned_masked'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_output('ndvi_base_aligned_masked'), + colormap=ndvi_colorramp + ), + # reclassified baseline lulc to ndvi values based on means + RasterPlotConfig( + raster_path=file_registry['ndvi_alt_aligned_masked'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_output('ndvi_alt_aligned_masked'), + colormap=ndvi_colorramp + ) + ]) + input_raster_config_list.append( + RasterPlotConfig( + raster_path=args_dict['population_raster'], + datatype=RasterDatatype.continuous, + spec=model_spec.get_input('population_raster') + ) + ) + + datatype = infer_continuous_or_divergent(file_registry['preventable_cases']) + is_continuous = datatype == RasterDatatype.continuous + common_kwargs = { + "datatype": datatype, + "colormap": 'Purples' if is_continuous else None, + "transform": RasterTransform.linear if is_continuous else RasterTransform.log, + "special_value_range": NEAR_ZERO_RANGE if is_continuous else None, + "special_value_color": '#FEE0B6' if is_continuous else None + } + output_raster_config_list = [ + SpecialRangeRasterPlotConfig( + raster_path=file_registry['preventable_cases'], + spec=model_spec.get_output('preventable_cases'), + **common_kwargs + ) + ] + if args_dict['health_cost_rate']: + output_raster_config_list.append( + SpecialRangeRasterPlotConfig( + raster_path=file_registry['preventable_cost'], + spec=model_spec.get_output('preventable_cost'), + **common_kwargs + ) + ) + + return (input_raster_config_list, + intermediate_raster_config_lists, + output_raster_config_list) + + +def _get_intermediate_output_headings(args_dict: dict) -> list[str]: + """Get headings for Intermediate Outputs sections of the report. + + Args: + args_dict (dict): the arguments passed to the model's ``execute`` + function. + + Returns: + A list containing exactly two strings or exactly three strings. + If the model was run with ``scenario==ndvi``, the report will display + a section with delta ndvi map and its inputs, and a second section + with baseline prevalence and cases. If the model was run with + ``scenario==lulc``, the report will also show the reclassified + LULC-to-NDVI rasters as intermediate outputs. + """ + intermediate_captions = [ + gettext('Difference in NDVI between Alternate and Baseline'), + gettext('Baseline Prevalence & Cases') + ] + if args_dict['scenario'] == 'lulc': + intermediate_captions.insert(0, + gettext('Reclassified Baseline & Alternate LULC (to NDVI)')) + + return intermediate_captions + + +def _generate_agg_results_table(file_registry: dict) -> str: + table_path = file_registry['preventable_cases_cost_sum_table'] + full_table_df = pandas.read_csv(table_path) + total_cases = list(full_table_df['total_cases'])[-1] + table_df = pandas.DataFrame({'Total Preventable Cases': [total_cases]}) + if file_registry.get('preventable_cost'): + total_cost = list(full_table_df['total_cost'])[-1] + table_df['Total Preventable Cost'] = [total_cost] + + return table_df.to_html(index=False) + + +def _create_aggregate_map( + geodataframe, + extent_feature, + xy_ratio, + attribute, + title): + """Create a choropleth map for a given attribute and return Vega JSON.""" + + # if the attribute has any negative values, use a divergent color + # scale; otherwise, use a continuous color scale + if (geodataframe[attribute] < 0).any(): + scale = altair.Scale(scheme='purpleorange', reverse=True, domainMid=0) + else: + scale = altair.Scale(scheme='purples') + + chart = altair.Chart(geodataframe).mark_geoshape( + clip=True, + stroke="white", + strokeWidth=0.5 + ).project( + type='identity', + reflectY=True, + fit=extent_feature + ).encode( + color=altair.Color( + f'{attribute}:Q', + scale=scale, + legend=altair.Legend(title=attribute) + ), + tooltip=[ + altair.Tooltip(f'{attribute}:Q', title=attribute, format=',.2f') + ] + ).properties( + width=MAP_WIDTH, + height=MAP_WIDTH / xy_ratio, + title=title + ).configure_legend(**vector_utils.LEGEND_CONFIG) + + return chart.to_json() + + +def report(file_registry: dict, args_dict: dict, model_spec: ModelSpec, + target_html_filepath: str): + """Generate an HTML summary of model results. + + Args: + file_registry (dict): The ``natcap.invest.FileRegistry.registry`` + that was returned by the model's ``execute`` method. + args_dict (dict): The arguments that were passed to the model's + ``execute`` method. + model_spec (natcap.invest.spec.ModelSpec): the model's ``MODEL_SPEC``. + target_html_filepath (str): path to an HTML file to be generated by + this function. + + Returns: + ``None`` + """ + + input_raster_config_list, \ + intermediate_raster_config_lists, \ + output_raster_config_list = _get_conditional_raster_plot_tuples( + model_spec, args_dict, file_registry) + + inputs_img_src = raster_utils.plot_and_base64_encode_rasters( + input_raster_config_list) + input_raster_caption = raster_utils.caption_raster_list( + input_raster_config_list) + + outputs_img_src = raster_utils.plot_and_base64_encode_rasters( + output_raster_config_list) + output_raster_caption = raster_utils.caption_raster_list( + output_raster_config_list) + + # There can be multiple sections for intermediate rasters + intermediate_img_srcs = [raster_utils.plot_and_base64_encode_rasters( + config_list) for config_list in intermediate_raster_config_lists] + intermediate_raster_captions = [raster_utils.caption_raster_list( + config_list) for config_list in intermediate_raster_config_lists] + + intermediate_headings = _get_intermediate_output_headings(args_dict) + + intermediate_raster_sections = [ + {'heading': heading, 'img_src': img_src, 'caption': caption} + for (heading, img_src, caption) + in zip(intermediate_headings, + intermediate_img_srcs, + intermediate_raster_captions) + ] + + input_raster_stats_table = raster_utils.raster_inputs_summary( + args_dict, model_spec).to_html(na_rep='') + + output_raster_stats_table = raster_utils.raster_workspace_summary( + file_registry).to_html(na_rep='') + + agg_results_table = _generate_agg_results_table(file_registry) + + # Vector maps + aggregate_vector_path = file_registry['preventable_cases_cost_sum_vector'] + aggregate_gdf = geopandas.read_file(aggregate_vector_path) + extent_feature, xy_ratio = vector_utils.get_geojson_bbox(aggregate_gdf) + + cases_map_json = _create_aggregate_map( + aggregate_gdf, + extent_feature, + xy_ratio, + 'sum_cases', + gettext('Preventable Cases by AOI Feature') + ) + cases_map_caption = [ + model_spec.get_output('preventable_cases_cost_sum_vector') + .get_field('sum_cases').about + ] + + cost_map_json = None + cost_map_caption = None + + if 'sum_cost' in aggregate_gdf.columns: + nonnull_cost = aggregate_gdf['sum_cost'].notna().any() + nonzero_cost = (aggregate_gdf['sum_cost'].fillna(0) != 0).any() + if nonnull_cost and nonzero_cost: + cost_map_json = _create_aggregate_map( + aggregate_gdf, + extent_feature, + xy_ratio, + 'sum_cost', + gettext('Preventable Cost by AOI Feature') + ) + cost_map_caption = [ + model_spec.get_output('preventable_cases_cost_sum_vector') + .get_field('sum_cost').about + ] + + aggregate_map_source_list = [ + model_spec.get_output('preventable_cases_cost_sum_vector').path + ] + + with open(target_html_filepath, 'w', encoding='utf-8') as target_file: + target_file.write(TEMPLATE.render( + report_script=model_spec.reporter, + invest_version=__version__, + report_filepath=target_html_filepath, + model_id=model_spec.model_id, + model_name=model_spec.model_title, + model_description=model_spec.about, + userguide_page=model_spec.userguide, + timestamp=time.strftime('%Y-%m-%d %H:%M'), + args_dict=args_dict, + agg_results_table=agg_results_table, + inputs_img_src=inputs_img_src, + inputs_caption=input_raster_caption, + outputs_img_src=outputs_img_src, + outputs_caption=output_raster_caption, + intermediate_raster_sections=intermediate_raster_sections, + raster_group_caption=report_constants.RASTER_GROUP_CAPTION, + lulc_pre_caption=report_constants.LULC_PRE_CAPTION, + output_raster_stats_table=output_raster_stats_table, + input_raster_stats_table=input_raster_stats_table, + stats_table_note=report_constants.STATS_TABLE_NOTE, + cases_map_json=cases_map_json, + cost_map_json=cost_map_json, + cases_map_caption=cases_map_caption, + cost_map_caption=cost_map_caption, + aggregate_map_source_list=aggregate_map_source_list, + model_spec_outputs=model_spec.outputs, + )) + + LOGGER.info(f'Created {target_html_filepath}') diff --git a/src/natcap/invest/urban_mental_health/urban_mental_health.py b/src/natcap/invest/urban_mental_health/urban_mental_health.py index ae5cd140a2..dc59c3ef8a 100644 --- a/src/natcap/invest/urban_mental_health/urban_mental_health.py +++ b/src/natcap/invest/urban_mental_health/urban_mental_health.py @@ -34,7 +34,7 @@ model_id="urban_mental_health", model_title=gettext("Urban Mental Health"), userguide="", # TODO - add this model to UG - reporter="", # TODO - add "natcap.invest.urban_mental_health.reporter", + reporter="natcap.invest.urban_mental_health.reporter", about=_model_description, validate_spatial_overlap=True, different_projections_ok=True, @@ -95,8 +95,8 @@ "increase in NDVI. If the user has an effect size value " "as an odds ratio, see User's Guide." ), - units=None, #todo: check - expression="value > 0" + units=None, + expression="value > 0 and value <= 1" ), spec.VectorInput( id="baseline_prevalence_vector", @@ -178,12 +178,13 @@ name=gettext("baseline land use/land cover"), about=gettext( "Map of land use/land cover codes under current or baseline " - "conditions. If scenario is not 'lulc', this is used only for " - "water masking. If an LULC attribute table is used, " - "all values in this raster must have corresponding entries. " - "This raster should extend beyond the AOI by at least the " - "search radius distance." - ), + "conditions. This raster should extend beyond the AOI by at " + "least the search radius distance. If an LULC attribute table " + "is provided, all values in this raster must have " + "corresponding entries. When using the NDVI scenario, this " + "raster may be used to mask out excluded land cover types " + "(such as water) based on an accompanying LULC attribute " + "table."), data_type=int, units=None, required="scenario=='lulc' or lulc_attr_csv", @@ -249,7 +250,9 @@ spec.SingleBandRasterOutput( id="preventable_cases", path="output/preventable_cases.tif", - about=gettext("Preventable cases at pixel level."), + about=gettext( + 'Preventable cases at pixel level. A negative value ' + 'indicates "excess" or "additional" cases.'), data_type=float, units=u.count, ), @@ -258,7 +261,8 @@ path="output/preventable_cost.tif", about=gettext( "Preventable cost at pixel level. The currency unit will " - "be the same as that in the health cost rate input."), + "be the same as that in the health cost rate input. A " + "negative value indicates an additional cost."), data_type=float, units=u.currency, created_if="health_cost_rate" @@ -328,14 +332,26 @@ spec.SingleBandRasterOutput( id="baseline_cases", path="intermediate/baseline_cases.tif", - about=gettext("Baseline cases raster."), + about=gettext( + "Baseline number of cases of the mental health outcome " + "per pixel, calculated as the rasterized baseline " + "prevalence rate multiplied by the population in each " + "pixel. This raster represents the expected number of " + "cases under baseline conditions and is used to estimate" + "preventable cases under the alternate/counterfactual " + "scenario."), data_type=float, units=u.count ), spec.SingleBandRasterOutput( id="baseline_prevalence_raster", path="intermediate/baseline_prevalence.tif", - about=gettext("Baseline prevalence raster."), + about=gettext( + "Rasterized baseline prevalence (or incidence) rate of " + "the mental health outcome. Pixel values are taken from " + "the `risk_rate` field of the baseline prevalence vector, " + "so each pixel within a polygon is assigned that " + "polygon's baseline rate."), data_type=float, units=None ), @@ -343,7 +359,7 @@ id="delta_ndvi", path="intermediate/delta_ndvi.tif", about=gettext( - "Difference between baseline and alternate NDVI raster."), + "Difference between alternate and baseline NDVI raster."), data_type=float, units=None ), @@ -435,7 +451,7 @@ "Preprocessed alternate NDVI raster. If using LULC " "inputs, this raster is created by masking, aligning, and " "resampling the alternate LULC and mapping it to mean " - "NDVI (with excluded lucodes set to NODATA)." + "NDVI (with excluded lucodes set to NODATA). " "If using NDVI inputs, this is simply the masked aligned " "and resampled alternate NDVI raster."), data_type=float, @@ -457,7 +473,7 @@ "Preprocessed baseline NDVI raster. If using LULC inputs, " "this raster is created by masking, aligning, and " "resampling the baseline LULC and mapping it to mean " - "NDVI (with excluded lucodes set to NODATA)." + "NDVI (with excluded lucodes set to NODATA). " "If using NDVI inputs, this is simply the masked aligned " "and resampled baseline NDVI raster."), data_type=float, @@ -887,7 +903,7 @@ def execute(args): ) zonal_stats_inputs = [ - args['aoi_path'], + args['aoi_path'], args['scenario'], file_registry['preventable_cases_cost_sum_table'], file_registry['preventable_cases_cost_sum_vector'], file_registry['preventable_cases']] @@ -1390,13 +1406,16 @@ def _preventable_cost_op(preventable_cases, cost, nodata): def zonal_stats_preventable_cases_cost( - base_vector_path, target_stats_csv, target_aggregate_vector_path, - preventable_cases_raster, preventable_cost_raster=None): + base_vector_path, scenario, target_stats_csv, + target_aggregate_vector_path, preventable_cases_raster, + preventable_cost_raster=None): """Calculate zonal statistics for each polygon in the AOI and write results to a csv and vector file. Args: base_vector_path (string): path to the AOI vector. + scenario (string): either 'ndvi' or 'lulc', used to label + output gpkg layer target_stats_csv (string): path to csv file to store dictionary returned by zonal stats. target_aggregate_vector_path (string): path to vector to store zonal @@ -1419,6 +1438,21 @@ def zonal_stats_preventable_cases_cost( driver.CreateCopy(target_aggregate_vector_path, aoi_vector) aoi_vector = None + # rename gpkg layer + custom_layer_name = f"preventable_cases_{scenario}" + + target_aggregate_vector = gdal.OpenEx( + target_aggregate_vector_path, gdal.OF_UPDATE) + + layer = target_aggregate_vector.GetLayer() + old_layer_name = layer.GetName() + layer = None + + target_aggregate_vector.ExecuteSQL( + f'ALTER TABLE "{old_layer_name}" RENAME TO "{custom_layer_name}"' + ) + target_aggregate_vector = None + cases_sum_field = ogr.FieldDefn('sum_cases', ogr.OFTReal) cases_sum_field.SetWidth(24) cases_sum_field.SetPrecision(11) diff --git a/tests/test_urban_mental_health.py b/tests/test_urban_mental_health.py index 2e0806ed9c..d3606e4c3c 100644 --- a/tests/test_urban_mental_health.py +++ b/tests/test_urban_mental_health.py @@ -732,7 +732,7 @@ def test_zonal_stats_preventable_cases_cost(self): target_agg_vector_path = os.path.join(self.workspace_dir, "stats.gpkg") urban_mental_health.zonal_stats_preventable_cases_cost( - aoi_vector, target_stats_csv, target_agg_vector_path, + aoi_vector, 'ndvi', target_stats_csv, target_agg_vector_path, preventable_cases, preventable_cost) # Check CSV