Update GEE SOC scripts for FarmLab (#9)

* Deleted indicators folder and edited README.md

* Added clhs workflow code and guide

* docs: clean up and formatting for cLHS readme

* chore: add a generally helpful gitignore file

* refactor: move clhs code into subfolder

* Force UTC time, added depth to BD calculation, added variables for SOC, BD and Label

---------

Co-authored-by: Ned Horning <[email protected]>
Co-authored-by: Kyle Lawlor-Bagcal <[email protected]>

socMapping/StockSOC_ExtractPoints.ipynb

@@ -81,6 +81,20 @@
     "outPickle = \"\""
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2d444d71",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "### Define the attribute labels from the input tabular data for SOC, BD, and the point name ###\n",
+    "# The attribute labels are the same as the attribute names in the point location ESRI Shapefile\n",
+    "SOC = 'C%'  # Attribute name for soil carbon metric \n",
+    "BD = 'BD'   # Attribute name for bulk density\n",
+    "depth = 15    # Soil sample depth in cm"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 5,
@@ -269,7 +283,7 @@
     "# Create a list of dates for all images in the collection\n",
     "datesObject = sentinelCollection.aggregate_array(\"system:time_start\")\n",
     "dateList =  datesObject.getInfo()\n",
-    "dateList=[datetime.fromtimestamp(x/1000).strftime('%Y_%m_%d') for x in dateList]"
+    "dateList = [datetime.utcfromtimestamp(x / 1000).strftime('%Y_%m_%d') for x in dateList]"
    ]
   },
   {
@@ -394,7 +408,7 @@
     "    dictarr = getValues(extractedPoints)\n",
     "    points = gpd.GeoDataFrame(dictarr)\n",
     "    # Add the following variables to the collection of point data\n",
-    "    points['stock'] = points['BD'] * points['C%']\n",
+    "    points['stock'] = points[BD] * points[SOC] * depth\n",
     "    points['twi'] = gpd.GeoDataFrame(dictarrTWI)['first']\n",
     "    points['chili'] = gpd.GeoDataFrame(dictarrCHILI)['first']\n",
     "    \n",
@@ -2206,7 +2220,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.10.8"
   }
  },
  "nbformat": 4,

socMapping/StockSOC_ExtractPoints.py

@@ -66,6 +66,16 @@
 outPickle = ""
 
 
+# In[ ]:
+
+
+### Define the attribute labels from the input tabular data for SOC, BD, and the point name ###
+# The attribute labels are the same as the attribute names in the point location ESRI Shapefile
+SOC = 'C%'  # Attribute name for soil carbon metric 
+BD = 'BD'   # Attribute name for bulk density
+depth = 15    # Soil sample depth in cm
+
+
 # In[5]:
 
 
@@ -214,7 +224,7 @@ def getValues(fc):
 # Create a list of dates for all images in the collection
 datesObject = sentinelCollection.aggregate_array("system:time_start")
 dateList =  datesObject.getInfo()
-dateList=[datetime.fromtimestamp(x/1000).strftime('%Y_%m_%d') for x in dateList]
+dateList = [datetime.utcfromtimestamp(x / 1000).strftime('%Y_%m_%d') for x in dateList]
 
 
 # In[15]:
@@ -287,7 +297,7 @@ def getValues(fc):
     dictarr = getValues(extractedPoints)
     points = gpd.GeoDataFrame(dictarr)
     # Add the following variables to the collection of point data
-    points['stock'] = points['BD'] * points['C%']
+    points['stock'] = points[BD] * points[SOC] * depth
     points['twi'] = gpd.GeoDataFrame(dictarrTWI)['first']
     points['chili'] = gpd.GeoDataFrame(dictarrCHILI)['first']
 

socMapping/StockSOC_ProcessPoints.ipynb

@@ -859,7 +859,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.10.8"
   }
  },
  "nbformat": 4,

socMapping/StockSOC_ProcessPoints.py

@@ -141,8 +141,8 @@ def adjust_r2(r2, num_examples, num_features):
 # In[15]:
 
 
-# Iterate through the dictionary of tabular one date at a time to find the set of variables
-# that gievs the lowest R2 value
+# Iterate through the dictionary one date at a time to find the set of variables
+# that gives the highest R2 value
 for iteration, key in enumerate(pointsDFs):
     points = pointsDFs[key]
     if (points['B3'].isna().sum() / len(points.index) < 0.2): # If under 20% of the values are NA (cloud masked)
@@ -178,7 +178,7 @@ def adjust_r2(r2, num_examples, num_features):
         for i in efs.subsets_:
             efs.subsets_[i]['adjusted_avg_score'] = (
             adjust_r2(r2=efs.subsets_[i]['avg_score'],
-                  num_examples=x.shape[0]/10,
+                  num_examples=x.shape[0]/1.0,
                   num_features=len(efs.subsets_[i]['feature_idx']))
             )
         score = -99e10

socMapping/stockSOC_PredictImage.ipynb

@@ -521,8 +521,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "# Multiply the output image by 1000 to be able to convert to integer allowing larger areas to be downloaded\n",
-    "outputImage = predImage.multiply(1000).round().toInt16()"
+    "# Multiply the output image by 10 to be able to convert to integer allowing larger areas to be downloaded\n",
+    "outputImage = predImage.multiply(10).round().toInt16()"
    ]
   },
   {
@@ -564,7 +564,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.9.13"
+   "version": "3.10.8"
   }
  },
  "nbformat": 4,

socMapping/stockSOC_PredictImage.py

@@ -368,8 +368,8 @@ def getValues(fc):
 # In[29]:
 
 
-# Multiply the output image by 1000 to be able to convert to integer allowing larger areas to be downloaded
-outputImage = predImage.multiply(1000).round().toInt16()
+# Multiply the output image by 10 to be able to convert to integer allowing larger areas to be downloaded
+outputImage = predImage.multiply(10).round().toInt16()
 
 
 # In[30]: