do several changes to improve the export function

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,

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",

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))

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)