From a2df0d673edf6af37df4e785fdab9bc488f43dbd Mon Sep 17 00:00:00 2001 From: cesarLima1 Date: Fri, 20 Dec 2024 11:57:20 -0400 Subject: [PATCH 1/2] [TM-1467] Adjust GFW data areas with correction factor using ESRI:54009 --- app/Services/RunIndicatorAnalysisService.php | 49 ++++++++- resources/python/gfw-area-adjustment/app.py | 101 +++++++++++++++++++ 2 files changed, 149 insertions(+), 1 deletion(-) create mode 100644 resources/python/gfw-area-adjustment/app.py diff --git a/app/Services/RunIndicatorAnalysisService.php b/app/Services/RunIndicatorAnalysisService.php index 6822441b..a8a60125 100644 --- a/app/Services/RunIndicatorAnalysisService.php +++ b/app/Services/RunIndicatorAnalysisService.php @@ -10,6 +10,7 @@ use Illuminate\Support\Facades\DB; use Illuminate\Support\Facades\Http; use Illuminate\Support\Facades\Log; +use Symfony\Component\Process\Process; class RunIndicatorAnalysisService { @@ -133,13 +134,59 @@ public function getGeometry($polygonUuid) public function sendApiRequestIndicator($secret_key, $query_url, $query_sql, $geometry) { - return Http::withHeaders([ + $response = Http::withHeaders([ 'content-type' => 'application/json', 'x-api-key' => $secret_key, ])->post('https://data-api.globalforestwatch.org' . $query_url, [ 'sql' => $query_sql, 'geometry' => $geometry, ]); + + if ($response->successful()) { + $gfwDataFile = tempnam(sys_get_temp_dir(), 'gfw_') . '.json'; + $geometryFile = tempnam(sys_get_temp_dir(), 'geom_') . '.json'; + $outputFile = tempnam(sys_get_temp_dir(), 'output_') . '.json'; + + try { + file_put_contents($gfwDataFile, json_encode($response->json())); + file_put_contents($geometryFile, json_encode($geometry)); + + $process = new Process([ + 'python3', + base_path() . '/resources/python/gfw-area-adjustment/app.py', + $gfwDataFile, + $geometryFile, + $outputFile, + ]); + + $process->run(); + + if (!$process->isSuccessful()) { + Log::error('Area adjustment failed: ' . $process->getErrorOutput()); + return $response; + } + + $adjustedData = json_decode(file_get_contents($outputFile), true); + + return new \Illuminate\Http\Client\Response( + new \GuzzleHttp\Psr7\Response( + 200, + ['Content-Type' => 'application/json'], + json_encode($adjustedData) + ) + ); + + } catch (\Exception $e) { + Log::error('Error adjusting areas: ' . $e->getMessage()); + return $response; + } finally { + @unlink($gfwDataFile); + @unlink($geometryFile); + @unlink($outputFile); + } + } + + return $response; } public function processTreeCoverLossValue($data, $indicator) diff --git a/resources/python/gfw-area-adjustment/app.py b/resources/python/gfw-area-adjustment/app.py new file mode 100644 index 00000000..b586a60d --- /dev/null +++ b/resources/python/gfw-area-adjustment/app.py @@ -0,0 +1,101 @@ +import sys +import json +import geopandas as gpd +from shapely.geometry import shape, Polygon +import logging + +logging.basicConfig(level=logging.INFO) +logger = logging.getLogger(__name__) + +def calculate_correction_factor(reference_geometry): + """ + Calculate the correction factor by comparing WGS84 and ESRI:54009 areas + All areas are converted to hectares (1 hectare = 10000 square meters) + """ + try: + gdf_wgs84 = gpd.GeoDataFrame(geometry=[reference_geometry], crs="EPSG:4326") + + gdf_projected = gdf_wgs84.to_crs('ESRI:54009') + area_projected_ha = gdf_projected.area[0] / 10000 + + if area_projected_ha < 0.0001: + logger.warning("Area is very small, using correction factor of 1") + return 1.0 + + area_geodesic_ha = gdf_wgs84.geometry.to_crs('+proj=cea').area[0] / 10000 + + correction_factor = area_projected_ha / area_geodesic_ha if area_geodesic_ha > 0 else 1 + + if correction_factor > 10 or correction_factor < 0.1: + logger.warning(f"Extreme correction factor detected: {correction_factor}. Using 1.0 instead.") + return 1.0 + + return correction_factor + + except Exception as e: + logger.error(f"Error calculating correction factor: {str(e)}") + raise + +def adjust_gfw_data(gfw_data, reference_geometry): + """ + Adjust GFW area values using the correction factor from reference geometry + """ + try: + if isinstance(reference_geometry, str): + reference_geometry = json.loads(reference_geometry) + + if 'type' in reference_geometry and reference_geometry['type'] == 'Feature': + geometry = shape(reference_geometry['geometry']) + elif 'type' in reference_geometry and reference_geometry['type'] == 'FeatureCollection': + geometry = shape(reference_geometry['features'][0]['geometry']) + else: + geometry = shape(reference_geometry) + + correction_factor = calculate_correction_factor(geometry) + + if isinstance(gfw_data, str): + gfw_data = json.loads(gfw_data) + + adjusted_data = { + "data": [], + "status": gfw_data.get("status", "success") + } + + for entry in gfw_data.get("data", []): + adjusted_entry = entry.copy() + if entry["area__ha"] > 0.0001: + adjusted_entry["area__ha"] = round(entry["area__ha"] * correction_factor, 5) + adjusted_data["data"].append(adjusted_entry) + + return adjusted_data + + except Exception as e: + logger.error(f"Error adjusting GFW data: {str(e)}") + raise + +def main(): + try: + if len(sys.argv) != 4: + raise ValueError("Script requires GFW data, reference geometry, and output file paths as arguments") + + gfw_data_file = sys.argv[1] + reference_geometry_file = sys.argv[2] + output_file = sys.argv[3] + + with open(gfw_data_file, 'r') as f: + gfw_data = json.load(f) + + with open(reference_geometry_file, 'r') as f: + reference_geometry = json.load(f) + + result = adjust_gfw_data(gfw_data, reference_geometry) + + with open(output_file, 'w') as f: + json.dump(result, f) + + except Exception as e: + logger.error(f"Error processing data: {str(e)}") + sys.exit(1) + +if __name__ == "__main__": + main() \ No newline at end of file From 3c1fe6338ac792f68b0c26ab68a0f33f464b5e05 Mon Sep 17 00:00:00 2001 From: cesarLima1 Date: Fri, 20 Dec 2024 12:20:43 -0400 Subject: [PATCH 2/2] [TM-1467] Adjust GFW data areas with correction factor using ESRI:54009 --- app/Services/RunIndicatorAnalysisService.php | 20 +++++++++++--------- 1 file changed, 11 insertions(+), 9 deletions(-) diff --git a/app/Services/RunIndicatorAnalysisService.php b/app/Services/RunIndicatorAnalysisService.php index a8a60125..19608be1 100644 --- a/app/Services/RunIndicatorAnalysisService.php +++ b/app/Services/RunIndicatorAnalysisService.php @@ -141,16 +141,16 @@ public function sendApiRequestIndicator($secret_key, $query_url, $query_sql, $ge 'sql' => $query_sql, 'geometry' => $geometry, ]); - + if ($response->successful()) { $gfwDataFile = tempnam(sys_get_temp_dir(), 'gfw_') . '.json'; $geometryFile = tempnam(sys_get_temp_dir(), 'geom_') . '.json'; $outputFile = tempnam(sys_get_temp_dir(), 'output_') . '.json'; - + try { file_put_contents($gfwDataFile, json_encode($response->json())); file_put_contents($geometryFile, json_encode($geometry)); - + $process = new Process([ 'python3', base_path() . '/resources/python/gfw-area-adjustment/app.py', @@ -158,14 +158,15 @@ public function sendApiRequestIndicator($secret_key, $query_url, $query_sql, $ge $geometryFile, $outputFile, ]); - + $process->run(); - - if (!$process->isSuccessful()) { + + if (! $process->isSuccessful()) { Log::error('Area adjustment failed: ' . $process->getErrorOutput()); + return $response; } - + $adjustedData = json_decode(file_get_contents($outputFile), true); return new \Illuminate\Http\Client\Response( @@ -175,9 +176,10 @@ public function sendApiRequestIndicator($secret_key, $query_url, $query_sql, $ge json_encode($adjustedData) ) ); - + } catch (\Exception $e) { Log::error('Error adjusting areas: ' . $e->getMessage()); + return $response; } finally { @unlink($gfwDataFile); @@ -185,7 +187,7 @@ public function sendApiRequestIndicator($secret_key, $query_url, $query_sql, $ge @unlink($outputFile); } } - + return $response; }