From 63f3566c8a009d60516fbaee6ea1524589b460c6 Mon Sep 17 00:00:00 2001 From: George Silva Date: Thu, 19 Dec 2024 18:08:28 -0300 Subject: [PATCH] do several changes to improve the export function --- src/planscape/planning/models.py | 17 +++++- src/planscape/planning/services.py | 4 +- src/planscape/planning/tests/test_services.py | 61 +++++++++++-------- src/planscape/rscripts/forsys.R | 36 ++++++----- 4 files changed, 76 insertions(+), 42 deletions(-) diff --git a/src/planscape/planning/models.py b/src/planscape/planning/models.py index fe10bd665..dbff4b7a2 100644 --- a/src/planscape/planning/models.py +++ b/src/planscape/planning/models.py @@ -1,3 +1,4 @@ +import json import uuid from pathlib import Path from typing import Optional, Type @@ -234,6 +235,18 @@ def get_forsys_folder(self) -> Path: def get_stand_size(self) -> StandSizeChoices: return self.configuration.get("stand_size", {}) or StandSizeChoices.LARGE + def get_geojson_result(self): + features = [ + { + "type": "Feature", + "properties": {**project_area.data, "name": project_area.name}, + "id": project_area.id, + "geometry": json.loads(project_area.geometry.geojson), + } + for project_area in self.project_areas.all() + ] + return {"type": "FeatureCollection", "features": features} + objects = ScenarioManager() class Meta(TypedModelMeta): @@ -335,7 +348,9 @@ class ProjectArea( name = models.CharField(max_length=128, help_text="Name of the Project Area.") - data = models.JSONField(null=True, help_text="Project Area data from Forsys.") + data = models.JSONField( + null=True, help_text="Project Area data from Forsys.", encoder=DjangoJSONEncoder + ) geometry = models.MultiPolygonField( srid=settings.CRS_INTERNAL_REPRESENTATION, diff --git a/src/planscape/planning/services.py b/src/planscape/planning/services.py index 366b82bfe..a3f20e5a5 100644 --- a/src/planscape/planning/services.py +++ b/src/planscape/planning/services.py @@ -349,14 +349,14 @@ def export_to_shapefile(scenario: Scenario) -> Path: if scenario.results.status != ScenarioResultStatus.SUCCESS: raise ValueError("Cannot export a scenario if it's result failed.") - geojson = scenario.results.result + geojson = scenario.get_geojson_result() schema = get_schema(geojson) shapefile_folder = scenario.get_shapefile_folder() shapefile_file = f"{scenario.name}.shp" shapefile_path = shapefile_folder / shapefile_file if not shapefile_folder.exists(): shapefile_folder.mkdir(parents=True) - crs = from_epsg(4326) + crs = from_epsg(settings.CRS_INTERNAL_REPRESENTATION) with fiona.open( str(shapefile_path), "w", diff --git a/src/planscape/planning/tests/test_services.py b/src/planscape/planning/tests/test_services.py index 61f266b30..3a3e1454e 100644 --- a/src/planscape/planning/tests/test_services.py +++ b/src/planscape/planning/tests/test_services.py @@ -3,6 +3,7 @@ from django.contrib.gis.geos import GEOSGeometry, MultiPolygon from django.test import TestCase, TransactionTestCase import fiona +from planscape.tests.factories import UserFactory from fiona.crs import to_string from planning.services import ( export_to_shapefile, @@ -11,7 +12,13 @@ get_schema, validate_scenario_treatment_ratio, ) -from planning.models import PlanningArea, Scenario, ScenarioResult, ScenarioResultStatus +from planning.models import ( + PlanningArea, + ProjectArea, + Scenario, + ScenarioResult, + ScenarioResultStatus, +) from stands.models import Stand, StandSizeChoices @@ -201,39 +208,43 @@ def test_export_raises_value_error_running(self): export_to_shapefile(scenario) def test_export_creates_file(self): + user = UserFactory.create() unit_poly = GEOSGeometry( "MULTIPOLYGON (((0 0, 0 1, 1 1, 1 0, 0 0)))", srid=4269 ) planning = PlanningArea.objects.create( - name="foo", region_name="sierra-nevada", geometry=unit_poly - ) - scenario = Scenario.objects.create(planning_area=planning, name="s1") - result = ScenarioResult.objects.create( + name="foo", + region_name="sierra-nevada", + geometry=unit_poly, + user=user, + ) + scenario = Scenario.objects.create( + planning_area=planning, + name="s1", + user=user, + ) + data = { + "foo": "abc", + "bar": 1, + "baz": 1.2, + "now": str(datetime.now()), + "today": date.today(), + } + ProjectArea.objects.create( scenario=scenario, - status=ScenarioResultStatus.SUCCESS, - result={ - "type": "FeatureCollection", - "features": [ - { - "properties": { - "foo": "abc", - "bar": 1, - "baz": 1.2, - "now": str(datetime.now()), - "today": date.today(), - }, - "geometry": { - "type": "Polygon", - "coordinates": [[[0, 0], [0, 1], [1, 1], [1, 0], [0, 0]]], - }, - } - ], - }, + geometry=unit_poly, + data=data, + created_by=user, + name="foo", ) + ScenarioResult.objects.create( + scenario=scenario, status=ScenarioResultStatus.SUCCESS + ) + output = export_to_shapefile(scenario) self.assertIsNotNone(output) path = output / "s1.shp" with fiona.open(path, "r", "ESRI Shapefile") as source: self.assertEqual(1, len(source)) - self.assertEqual(to_string(source.crs), "EPSG:4326") + self.assertEqual(to_string(source.crs), "EPSG:4269") shutil.rmtree(str(output)) diff --git a/src/planscape/rscripts/forsys.R b/src/planscape/rscripts/forsys.R index 88ab9cbb0..64bd4e0cc 100644 --- a/src/planscape/rscripts/forsys.R +++ b/src/planscape/rscripts/forsys.R @@ -376,11 +376,7 @@ get_stand_metrics <- function( get_project_geometry <- function(connection, stand_ids) { query <- glue_sql("SELECT ST_AsGeoJSON( - ST_Transform( - ST_Union(geometry), - 4326), - 6, -- max precision - 0 -- output mode + ST_Union(geometry) ) FROM stands_stand WHERE id IN ({stand_ids*})", @@ -391,6 +387,19 @@ get_project_geometry <- function(connection, stand_ids) { return(fromJSON(result$st_asgeojson)) } +get_project_geometry_text <- function(connection, stand_ids) { + query <- glue_sql("SELECT + ST_AsText( + ST_Union(geometry) + ) as \"geometry\" + FROM stands_stand + WHERE id IN ({stand_ids*})", + stand_ids = stand_ids, + .con = connection + ) + result <- dbGetQuery(connection, query) + return(result$geometry) +} get_project_ids <- function(forsys_output) { return(unique(forsys_output$project_output$proj_id)) @@ -422,7 +431,8 @@ to_properties <- function( forsys_project_outputs, project_stand_count, stand_size, - new_column_for_postprocessing = FALSE) { + new_column_for_postprocessing = FALSE, + text_geometry) { scenario_cost_per_acre <- get_cost_per_acre(scenario) project_data <- forsys_project_outputs %>% filter(proj_id == project_id) %>% @@ -431,6 +441,7 @@ to_properties <- function( mutate(total_cost = ETrt_area_acres * scenario_cost_per_acre) %>% mutate(cost_per_acre = scenario_cost_per_acre) %>% mutate(pct_area = ETrt_area_acres / scenario$planning_area_acres) %>% + mutate(text_geometry = text_geometry) %>% rename_with(.fn = rename_col) # post process @@ -473,13 +484,15 @@ to_project_data <- function( stand_count <- nrow(project_stand_ids) project_stand_ids <- as.integer(project_stand_ids$stand_id) geometry <- get_project_geometry(connection, project_stand_ids) + text_geometry <- get_project_geometry_text(connection, project_stand_ids) properties <- to_properties( project_id, scenario, forsys_outputs$project_output, stand_count, stand_size, - new_column_for_postprocessing + new_column_for_postprocessing, + text_geometry ) return(list( type = "Feature", @@ -835,12 +848,7 @@ upsert_project_area <- function( {data}, -- we have to convert to 4269 because ST_GeomFromGeoJSON assumes 4326 ST_Multi( - ST_Transform( - ST_GeomFromGeoJSON( - {geometry} - ), - 4269 - ) + ST_SetSRID(ST_GeomFromText({geometry}), 4269) ) ) ON CONFLICT (scenario_id, name) DO UPDATE @@ -858,7 +866,7 @@ upsert_project_area <- function( scenario_id=scenario$id, name=area_name, data=toJSON(project$properties), - geometry=toJSON(project$geometry), + geometry=project$properties$text_geometry, .con = connection ) dbExecute(connection, query, immediate = TRUE)