fix: traceback errors-geo

jjmaynard · jjmaynard · commit 40261a0218ed · 2025-03-28T07:53:50.000-07:00
-fixed distance calculation errors associated with coordinate system transformation.
-added top depth input to rank function
diff --git a/Makefile b/Makefile
@@ -57,10 +57,10 @@ process_bulk_test_results_us:
 	python -m soil_id.tests.us.process_bulk_test_results $(RESULTS_FILE)
 
 generate_bulk_test_results_global:
-	python -m soil_id.tests.us.generate_bulk_test_results
+	python -m soil_id.tests.global.generate_bulk_test_results
 
 process_bulk_test_results_global:
-	python -m soil_id.tests.us.process_bulk_test_results $(RESULTS_FILE)
+	python -m soil_id.tests.global.process_bulk_test_results $(RESULTS_FILE)
 
 # Donwload Munsell CSV, SHX, SHP, SBX, SBN, PRJ, DBF
 download-soil-data:
diff --git a/soil_id/db.py b/soil_id/db.py
@@ -27,9 +27,6 @@
 # local libraries
 import soil_id.config
 
-from .utils import get_target_utm_srid
-
-
 def get_datastore_connection():
     """
     Establish a connection to the datastore using app configurations.
@@ -208,9 +205,11 @@ def get_hwsd2_profile_data(conn, hwsd2_mu_select):
         logging.error(f"Error querying PostgreSQL: {err}")
         return pd.DataFrame()
 
+
 def extract_hwsd2_data(lon, lat, buffer_dist, table_name):
     """
-    Fetches HWSD soil data from a PostGIS table within a given buffer around a point.
+    Fetches HWSD soil data from a PostGIS table within a given buffer around a point,
+    performing distance and intersection calculations directly on geographic coordinates.
     
     Parameters:
         lon (float): Longitude of the problem point.
@@ -221,35 +220,30 @@ def extract_hwsd2_data(lon, lat, buffer_dist, table_name):
     Returns:
         DataFrame: Merged data from hwsdv2 and hwsdv2_data.
     """
-    # Determine target UTM SRID as an integer (e.g., 32642)
-    target_srid = get_target_utm_srid(lat, lon)
-    
     # Use a single connection for both queries.
     with get_datastore_connection() as conn:
         # Compute the buffer polygon (in WKT) around the problem point.
+        # Here, we use the geography type to compute a buffer in meters,
+        # then cast it back to geometry in EPSG:4326.
         buffer_query = """
             WITH buffer AS (
               SELECT ST_AsText(
-                       ST_Transform(
-                         ST_Buffer(
-                           ST_Transform(
-                             ST_SetSRID(ST_Point(%s, %s), 4326),
-                             %s
-                           ),
-                           %s
-                         ),
-                         4326
-                       )
+                       ST_Buffer(
+                         ST_SetSRID(ST_Point(%s, %s), 4326)::geography,
+                         %s
+                       )::geometry
                      ) AS wkt
             )
             SELECT wkt FROM buffer;
         """
         with conn.cursor() as cur:
-            cur.execute(buffer_query, (lon, lat, target_srid, buffer_dist))
+            cur.execute(buffer_query, (lon, lat, buffer_dist))
             buffer_wkt = cur.fetchone()[0]
             print("Buffer WKT:", buffer_wkt)
         
         # Build the main query that uses the computed buffer.
+        # Distance is computed by casting geometries to geography,
+        # which returns the geodesic distance in meters.
         main_query = f"""
             WITH valid_geom AS (
               SELECT
@@ -261,19 +255,19 @@ def extract_hwsd2_data(lon, lat, buffer_dist, table_name):
             )
             SELECT
               hwsd2,
-              ST_AsEWKB(ST_Transform(geom, {target_srid})) AS geom,
+              ST_AsEWKB(geom) AS geom,
               ST_Distance(
-                ST_Transform(geom, {target_srid}),
-                ST_Transform(ST_SetSRID(ST_Point({lon}, {lat}), 4326), {target_srid})
+                geom::geography,
+                ST_SetSRID(ST_Point({lon}, {lat}), 4326)::geography
               ) AS distance,
               ST_Intersects(
-                ST_Transform(geom, {target_srid}),
-                ST_Transform(ST_SetSRID(ST_Point({lon}, {lat}), 4326), {target_srid})
+                geom,
+                ST_SetSRID(ST_Point({lon}, {lat}), 4326)
               ) AS pt_intersect
             FROM valid_geom
             WHERE ST_Intersects(
-                ST_Transform(geom, {target_srid}),
-                ST_Transform(ST_SetSRID(ST_Point({lon}, {lat}), 4326), {target_srid})
+                geom,
+                ST_SetSRID(ST_Point({lon}, {lat}), 4326)
             );
         """
         
diff --git a/soil_id/global_soil.py b/soil_id/global_soil.py
@@ -79,12 +79,15 @@ class SoilListOutputData:
 ##################################################################################################
 def list_soils_global(lon, lat):
     # Extract HWSD2 Data
-    hwsd2_data = extract_hwsd2_data(
-        lon,
-        lat,
-        table_name="hwsdv2",
-        buffer_dist=10000,
-    )
+    try:
+        hwsd2_data = extract_hwsd2_data(
+            lon,
+            lat,
+            table_name="hwsdv2",
+            buffer_dist=10000,
+        )
+    except KeyError:
+        return("Soil map information is not available at this location")
 
     # Component Data
     mucompdata_pd = hwsd2_data[["hwsd2", "fao90_name", "distance", "share", "compid"]]
@@ -556,7 +559,8 @@ def rank_soils_global(
     lat,
     list_output_data: SoilListOutputData,
     soilHorizon,
-    horizonDepth,
+    topDepth,
+    bottomDepth,
     rfvDepth,
     lab_Color,
     bedrock,
@@ -568,7 +572,8 @@ def rank_soils_global(
     soil_df = pd.DataFrame(
         {
             "soilHorizon": soilHorizon,
-            "horizonDepth": horizonDepth,
+            "top": topDepth,
+            "bottom": bottomDepth,
             "rfvDepth": rfvDepth,
             "lab_Color": lab_Color,
         }
@@ -577,15 +582,9 @@ def rank_soils_global(
     # Drop rows where all values are NaN
     soil_df.dropna(how="all", inplace=True)
 
-    # Set the bottom of each horizon
-    soil_df["bottom"] = soil_df["horizonDepth"]
-
     # Replace NaNs with None for consistency
     # soil_df.fillna(value=None, inplace=True)
 
-    # Calculate the top depth for each horizon
-    soil_df["top"] = [0] + soil_df["horizonDepth"].iloc[:-1].tolist()
-
     # Adjust the bottom depth based on bedrock depth
     if bedrock is not None:
         if bedrock is not soil_df.empty and soil_df["bottom"].iloc[-1] > bedrock:
@@ -596,9 +595,6 @@ def rank_soils_global(
             # Set the bottom depth of the last row to the bedrock depth
             soil_df.at[last_valid_index, "bottom"] = bedrock
 
-    # Drop the original horizonDepth column
-    soil_df.drop(columns=["horizonDepth"], inplace=True)
-
     # Filter out rows without valid horizon data
     relevant_columns = ["soilHorizon", "rfvDepth", "lab_Color"]
     soil_df_slice = soil_df.dropna(how="all", subset=relevant_columns)
@@ -623,8 +619,8 @@ def rank_soils_global(
         # Convert soil properties to lists
         soilHorizon = soil_df.soilHorizon.tolist()
         rfvDepth = soil_df.rfvDepth.tolist()
-        horizonDepthB = [int(x) for x in soil_df.bottom.tolist()]
-        horizonDepthT = [int(x) for x in soil_df.top.tolist()]
+        bottom = [int(x) for x in soil_df.bottom.tolist()]
+        top = [int(x) for x in soil_df.top.tolist()]
         lab_Color = soil_df.lab_Color
 
         # Generate user specified percent clay, sand, and rfv distributions
@@ -635,17 +631,17 @@ def rank_soils_global(
         p_sandpct_intpl = [
             spt[i]
             for i in range(len(soilHorizon))
-            for _ in range(horizonDepthT[i], horizonDepthB[i])
+            for _ in range(top[i], bottom[i])
         ]
         p_claypct_intpl = [
             cpt[i]
             for i in range(len(soilHorizon))
-            for _ in range(horizonDepthT[i], horizonDepthB[i])
+            for _ in range(top[i], bottom[i])
         ]
         p_cfg_intpl = [
             p_cfg[i]
             for i in range(len(soilHorizon))
-            for _ in range(horizonDepthT[i], horizonDepthB[i])
+            for _ in range(top[i], bottom[i])
         ]
 
         # Length of interpolated texture and RF depth
@@ -657,7 +653,7 @@ def rank_soils_global(
         if not lab_Color.isnull().all():
             lab_Color = [[np.nan, np.nan, np.nan] if x is None else x for x in lab_Color]
             lab_Color = pd.DataFrame(lab_Color)
-            pedon_LAB = pedon_color(lab_Color, horizonDepth)
+            pedon_LAB = pedon_color(lab_Color, top, bottom)
 
             if not np.isnan(pedon_LAB).all():
                 refs = {
@@ -840,8 +836,8 @@ def rank_soils_global(
         dis_max = max(map(np.nanmax, dis_mat_list))
 
         # Apply depth weight
-        depth_weight = np.concatenate([np.repeat(0.2, 20), np.repeat(1.0, 80)])
-        depth_weight = depth_weight[: len(soil_matrix)]
+        depth_weight = np.concatenate([np.repeat(0.2, 20), np.repeat(1.0, 180)])
+        depth_weight = depth_weight[soil_matrix.index]
 
         # Infill NaN data
         for idx, dis_mat in enumerate(dis_mat_list):
diff --git a/soil_id/tests/global/generate_bulk_test_results.py b/soil_id/tests/global/generate_bulk_test_results.py
@@ -60,7 +60,8 @@
                     lat=lat,
                     lon=lon,
                     list_output_data=list_result,
-                    horizonDepth=pedon["BOTDEP"].values.tolist(),
+                    topDepth=pedon["TOPDEP"].values.tolist(),
+                    bottomDepth=pedon["BOTDEP"].values.tolist(),
                     soilHorizon=pedon["textClass"].values.tolist(),
                     rfvDepth=pedon["RFV"].values.tolist(),
                     lab_Color=pedon[["L", "A", "B"]].values.tolist(),
diff --git a/soil_id/tests/global/test_global.py b/soil_id/tests/global/test_global.py
@@ -38,7 +38,8 @@ def test_soil_location():
             item["lat"],
             list_output_data=list_soils_result,
             soilHorizon=["Loam"],
-            horizonDepth=[15],
+            topDepth=[15],
+            bottomDepth=[45],
             rfvDepth=[20],
             lab_Color=[[41.23035939, 3.623018224, 13.27654356]],
             bedrock=None,
diff --git a/soil_id/utils.py b/soil_id/utils.py
