[TM-1467] polygon area discrepancy

app/Console/Commands/RecalculatePolygonAreas.php

@@ -0,0 +1,94 @@
+<?php
+
+namespace App\Console\Commands;
+
+use App\Models\V2\PolygonGeometry;
+use App\Models\V2\Sites\SitePolygon;
+use App\Services\AreaCalculationService;
+use Illuminate\Console\Command;
+use Illuminate\Support\Facades\DB;
+
+class RecalculatePolygonAreas extends Command
+{
+    /**
+     * The name and signature of the console command.
+     *
+     * @var string
+     */
+    protected $signature = 'polygons:recalculate-areas 
+                            {--batch-size=100 : Number of records to process in each batch}
+                            {--only-active : Process only active polygons}';
+
+    /**
+     * The console command description.
+     *
+     * @var string
+     */
+    protected $description = 'Recalculate areas for all site polygons';
+
+    /**
+     * Execute the console command.
+     */
+    public function handle(AreaCalculationService $areaService)
+    {
+        DB::beginTransaction();
+
+        try {
+            $query = SitePolygon::query();
+
+            if ($this->option('only-active')) {
+                $query->active();
+            }
+
+            $sitePolygons = $query->cursor();
+
+            $processedCount = 0;
+            $errorCount = 0;
+
+            $this->info('Starting polygon area recalculation...');
+            $progressBar = $this->output->createProgressBar();
+            $progressBar->start();
+
+            foreach ($sitePolygons as $sitePolygon) {
+                try {
+                    $polygonGeometry = PolygonGeometry::where('uuid', $sitePolygon->poly_id)
+                    ->select('uuid', DB::raw('ST_AsGeoJSON(geom) AS geojsonGeometry'))
+                    ->first();
+                    if (! $polygonGeometry) {
+                        $this->error("No geometry found for poly_id: {$sitePolygon->poly_id}");
+                        $errorCount++;
+
+                        continue;
+                    }
+                    $geometry = json_decode($polygonGeometry->geojsonGeometry, true);
+
+                    $calculatedArea = $areaService->getArea($geometry);
+
+                    $sitePolygon->calc_area = $calculatedArea;
+                    $sitePolygon->save();
+
+                    $processedCount++;
+                    $progressBar->advance();
+                } catch (\Exception $e) {
+                    $this->error("Error processing polygon {$sitePolygon->id}: " . $e->getMessage());
+                    $errorCount++;
+                }
+            }
+
+            DB::commit();
+
+            $progressBar->finish();
+            $this->info("\n\nRecalculation complete!");
+            $this->info("Processed: {$processedCount} polygons");
+            $this->info("Errors: {$errorCount}");
+
+        } catch (\Exception $e) {
+            DB::rollBack();
+            $this->error('Recalculation failed: ' . $e->getMessage());
+
+            return self::FAILURE;
+        }
+
+        return self::SUCCESS;
+    }
+}

app/Helpers/PolygonGeometryHelper.php

@@ -5,7 +5,7 @@
 use App\Models\V2\Projects\Project;
 use App\Models\V2\Sites\Site;
 use App\Models\V2\Sites\SitePolygon;
-use Illuminate\Support\Facades\DB;
+use App\Services\AreaCalculationService;
 use Illuminate\Support\Facades\Log;
 
 class PolygonGeometryHelper
@@ -16,12 +16,8 @@ public static function updateEstAreainSitePolygon($polygonGeometry, $geometry)
             $sitePolygon = SitePolygon::where('poly_id', $polygonGeometry->uuid)->first();
 
             if ($sitePolygon) {
-                $geojson = json_encode($geometry);
-                $areaSqDegrees = DB::selectOne("SELECT ST_Area(ST_GeomFromGeoJSON('$geojson')) AS area")->area;
-                $latitude = DB::selectOne("SELECT ST_Y(ST_Centroid(ST_GeomFromGeoJSON('$geojson'))) AS latitude")->latitude;
-                $unitLatitude = 111320;
-                $areaSqMeters = $areaSqDegrees * pow($unitLatitude * cos(deg2rad($latitude)), 2);
-                $areaHectares = $areaSqMeters / 10000;
+                $areaCalculationService = app(AreaCalculationService::class);
+                $areaHectares = $areaCalculationService->getArea((array) $geometry->geometry);
 
                 $sitePolygon->calc_area = $areaHectares;
                 $sitePolygon->save();

app/Http/Controllers/V2/Terrafund/TerrafundEditGeometryController.php

@@ -9,6 +9,7 @@
 use App\Models\V2\Projects\ProjectPolygon;
 use App\Models\V2\Sites\SitePolygon;
 use App\Models\V2\User;
+use App\Services\AreaCalculationService;
 use App\Services\PolygonService;
 use App\Services\SiteService;
 use Illuminate\Http\Request;
@@ -329,10 +330,12 @@ public function createSitePolygon(string $uuid, string $siteUuid, Request $reque
             if (! $polygonGeometry) {
                 return response()->json(['message' => 'No polygon geometry found for the given UUID.'], 404);
             }
-            $areaSqDegrees = DB::selectOne('SELECT ST_Area(geom) AS area FROM polygon_geometry WHERE uuid = :uuid', ['uuid' => $uuid])->area;
-            $latitude = DB::selectOne('SELECT ST_Y(ST_Centroid(geom)) AS latitude FROM polygon_geometry WHERE uuid = :uuid', ['uuid' => $uuid])->latitude;
-            $areaSqMeters = $areaSqDegrees * pow(111320 * cos(deg2rad($latitude)), 2);
-            $areaHectares = $areaSqMeters / 10000;
+            $polygonGeom = PolygonGeometry::where('uuid', $uuid)
+            ->select('uuid', DB::raw('ST_AsGeoJSON(geom) AS geojsonGeometry'))
+            ->first();
+            $geometry = json_decode($polygonGeom->geojsonGeometry, true);
+            $areaCalculationService = app(AreaCalculationService::class);
+            $areaHectares = $areaCalculationService->getArea($geometry);
             $sitePolygon = new SitePolygon([
                 'poly_name' => $validatedData['poly_name'],
                 'plantstart' => $validatedData['plantstart'],

app/Services/AreaCalculationService.php

@@ -0,0 +1,85 @@
+<?php
+
+namespace App\Services;
+
+use Illuminate\Support\Facades\DB;
+use Illuminate\Support\Facades\Log;
+use Symfony\Component\Process\Process;
+
+class AreaCalculationService
+{
+    protected function calculateArea(array $geometry): float
+    {
+        $geojson = json_encode([
+            'type' => 'Feature',
+            'geometry' => $geometry,
+            'crs' => ['type' => 'name', 'properties' => ['name' => 'EPSG:4326']],
+        ]);
+
+        $inputGeojson = tempnam(sys_get_temp_dir(), 'input_') . '.geojson';
+        $outputGeojson = tempnam(sys_get_temp_dir(), 'output_') . '.geojson';
+
+        try {
+            file_put_contents($inputGeojson, $geojson);
+
+            $process = new Process([
+                'python3',
+                base_path() . '/resources/python/polygon-area/app.py',
+                $inputGeojson,
+                $outputGeojson,
+            ]);
+
+            $process->run();
+
+            if (! $process->isSuccessful()) {
+                Log::error('Area calculation failed: ' . $process->getErrorOutput());
+
+                throw new \RuntimeException('Area calculation failed: ' . $process->getErrorOutput());
+            }
+
+            $result = json_decode(file_get_contents($outputGeojson), true);
+
+            return $result['area_hectares'];
+
+        } catch (\Exception $e) {
+            Log::error('Error calculating area: ' . $e->getMessage());
+
+            throw $e;
+        } finally {
+            @unlink($inputGeojson);
+            @unlink($outputGeojson);
+        }
+    }
+
+    public function getGeomAndArea(array $geometry): array
+    {
+        $geojson = json_encode([
+            'type' => 'Feature',
+            'geometry' => $geometry,
+            'crs' => ['type' => 'name', 'properties' => ['name' => 'EPSG:4326']],
+        ]);
+
+        $geom = DB::raw("ST_GeomFromGeoJSON('$geojson')");
+        $areaHectares = $this->calculateArea($geometry);
+
+        return ['geom' => $geom, 'area' => $areaHectares];
+    }
+
+    public function getArea(array $geometry): float
+    {
+        if ($geometry['type'] === 'MultiPolygon') {
+            $totalArea = 0;
+            foreach ($geometry['coordinates'] as $polygon) {
+                $polygonGeometry = [
+                    'type' => 'Polygon',
+                    'coordinates' => $polygon,
+                ];
+                $totalArea += $this->calculateArea($polygonGeometry);
+            }
+
+            return $totalArea;
+        }
+
+        return $this->calculateArea($geometry);
+    }
+}

app/Services/PolygonService.php

@@ -284,20 +284,9 @@ protected function getGeom(array $geometry)
 
     protected function getGeomAndArea(array $geometry): array
     {
-        // Convert geometry to GeoJSON string
-        $geojson = json_encode(['type' => 'Feature', 'geometry' => $geometry, 'crs' => ['type' => 'name', 'properties' => ['name' => 'EPSG:4326']]]);
-
-        // Get GeoJSON data in the database
-        $geom = DB::raw("ST_GeomFromGeoJSON('$geojson')");
-        $areaSqDegrees = DB::selectOne("SELECT ST_Area(ST_GeomFromGeoJSON('$geojson')) AS area")->area;
-        $latitude = DB::selectOne("SELECT ST_Y(ST_Centroid(ST_GeomFromGeoJSON('$geojson'))) AS latitude")->latitude;
-        // 111320 is the length of one degree of latitude in meters at the equator
-        $unitLatitude = 111320;
-        $areaSqMeters = $areaSqDegrees * pow($unitLatitude * cos(deg2rad($latitude)), 2);
-
-        $areaHectares = $areaSqMeters / 10000;
+        $areaCalculationService = app(AreaCalculationService::class);
 
-        return ['geom' => $geom, 'area' => $areaHectares];
+        return $areaCalculationService->getGeomAndArea($geometry);
     }
 
     protected function insertSinglePolygon(array $geometry): array

docker/php.Dockerfile

@@ -1,7 +1,11 @@
 ## PHP
 FROM php:8.2-apache AS php
-RUN apt-get update
-RUN apt-get install -y \
+
+# Set GDAL version
+ENV GDAL_VERSION=3.4.3
+
+# Install basic dependencies
+RUN apt-get update && apt-get install -y \
     libxml2-dev \
     libonig-dev \
     libpng-dev \
@@ -10,10 +14,36 @@ RUN apt-get install -y \
     libmagickwand-dev \
     mariadb-client \
     libzip-dev \
-    gdal-bin \
-    libgdal-dev \
     python3.11-venv \
-    exiftool 
+    python3.11-dev \
+    exiftool \
+    build-essential \
+    wget \
+    cmake \
+    sqlite3 \
+    libsqlite3-dev \
+    libspatialite-dev \
+    libpq-dev \
+    libcurl4-gnutls-dev \
+    libproj-dev \
+    libgeos-dev \
+    && rm -rf /var/lib/apt/lists/*
+
+# Install GDAL 3.4.3 from source
+RUN wget https://github.com/OSGeo/gdal/releases/download/v${GDAL_VERSION}/gdal-${GDAL_VERSION}.tar.gz \
+    && tar xzf gdal-${GDAL_VERSION}.tar.gz \
+    && cd gdal-${GDAL_VERSION} \
+    && ./configure \
+    && make -j$(nproc) \
+    && make install \
+    && ldconfig \
+    && cd .. \
+    && rm -rf gdal-${GDAL_VERSION} gdal-${GDAL_VERSION}.tar.gz
+
+# Set GDAL environment variables
+ENV CPLUS_INCLUDE_PATH=/usr/include/gdal
+ENV C_INCLUDE_PATH=/usr/include/gdal
+ENV LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
 
 RUN docker-php-ext-configure gd --with-freetype --with-jpeg
 RUN docker-php-ext-install \
@@ -40,12 +70,16 @@ RUN a2enmod rewrite
 COPY docker/000-default.conf /etc/apache2/sites-available/000-default.conf
 COPY docker/php.ini /usr/local/etc/php/php.ini
 
-## Python
-RUN python3 -m venv /opt/python
-COPY resources/python/polygon-voronoi/requirements.txt /root/voronoi-requirements.txt
+# Python virtual environment setup
+RUN python3.11 -m venv /opt/python
 ENV PATH="/opt/python/bin:${PATH}"
+
+# Install Python dependencies in the correct order
+COPY resources/python/polygon-voronoi/requirements.txt /root/voronoi-requirements.txt
+RUN pip3 install --upgrade pip wheel setuptools
 RUN pip3 install -r /root/voronoi-requirements.txt
+
 RUN chmod -R a+rx /opt/python
 USER www-data
 ENV PATH="/opt/python/bin:${PATH}"
-USER root
+USER root

resources/python/polygon-area/app.py

@@ -0,0 +1,60 @@
+import sys
+import json
+import geopandas as gpd
+from shapely.geometry import shape
+import logging
+
+logging.basicConfig(level=logging.INFO)
+logger = logging.getLogger(__name__)
+
+def calculate_area(geometry):
+    """
+    Calculate area in hectares for a given geometry
+    """
+    try:
+        gdf = gpd.GeoDataFrame(geometry=[geometry], crs="EPSG:4326")
+        gdf_projected = gdf.to_crs('ESRI:54009')
+
+        area_hectares = gdf_projected.geometry.area[0] / 10000
+
+        return area_hectares
+    except Exception as e:
+        logger.error(f"Error calculating area: {str(e)}")
+        raise
+
+def main():
+    try:
+        if len(sys.argv) != 3:
+            raise ValueError("Script requires input and output file paths as arguments")
+
+        input_file = sys.argv[1]
+        output_file = sys.argv[2]
+
+        with open(input_file, 'r') as f:
+            geojson_data = json.load(f)
+
+        if 'type' in geojson_data and geojson_data['type'] == 'Feature':
+            geometry = shape(geojson_data['geometry'])
+        elif 'type' in geojson_data and geojson_data['type'] == 'FeatureCollection':
+            geometry = shape(geojson_data['features'][0]['geometry'])
+        else:
+            geometry = shape(geojson_data)
+
+        area = calculate_area(geometry)
+
+        result = {
+            'area_hectares': area,
+            'original_geometry': geojson_data
+        }
+
+        with open(output_file, 'w') as f:
+            json.dump(result, f)
+
+        print(area)
+
+    except Exception as e:
+        logger.error(f"Error processing geometry: {str(e)}")
+        sys.exit(1)
+
+if __name__ == "__main__":
+    main()