Merge pull request #2263 from WikiWatershed/arr/estimate-soilpals+bfigrd

Estimate SedPhos And Recession Coefficient From Rasters

src/mmw/apps/core/permissions.py

@@ -5,8 +5,6 @@
 from rest_framework.permissions import BasePermission
 from rest_framework.authentication import TokenAuthentication
 
-from django.conf import settings
-
 
 class IsTokenAuthenticatedOrNotSwagger(BasePermission):
     """

src/mmw/apps/modeling/mapshed/calcs.py

@@ -509,36 +509,6 @@ def cni_avg(nlcds):
     ]
 
 
-def sediment_phosphorus(nt_count):
-    """
-    Given a dictionary mapping pairs of NLCD Codes and Soil Textures to counts
-    of cells, returns the average concentration of phosphorus in the soil by
-    looking up the value in the SOILP table. NLCD Codes 81 and 82 correspond to
-    agricultural values, and the rest to non-agricultural ones.
-
-    Original at Class1.vb@1.3.0:8975-8988
-    """
-
-    ag_textures = {}
-    nag_textures = {}
-    total = sum(nt_count.values())
-
-    for (n, t), count in nt_count.iteritems():
-        if n in AG_NLCD_CODES:
-            ag_textures[t] = count + ag_textures.get(t, 0)
-        else:
-            nag_textures[t] = count + nag_textures.get(t, 0)
-
-    ag_sedp = sum(count * settings.SOILP[t][0]
-                  for t, count in ag_textures.iteritems())
-    nag_sedp = sum(count * settings.SOILP[t][1]
-                   for t, count in nag_textures.iteritems())
-
-    sedp = float(ag_sedp + nag_sedp) / total
-
-    return sedp * 1.6
-
-
 def groundwater_nitrogen_conc(gwn_dict):
     """
     Calculate GwN and GwP from a dictionary of NLCD land use keys

src/mmw/apps/modeling/mapshed/tasks.py

@@ -4,7 +4,6 @@
 from __future__ import absolute_import
 
 from celery import shared_task
-
 from django.conf import settings
 from django.contrib.gis.geos import GEOSGeometry
 
@@ -25,7 +24,6 @@
                                          point_source_discharge,
                                          weather_data,
                                          curve_number,
-                                         sediment_phosphorus,
                                          phosphorus_conc,
                                          groundwater_nitrogen_conc,
                                          sediment_delivery_ratio,
@@ -106,7 +104,7 @@ def collect_data(geop_result, geojson):
     z['n46f'] = geop_result['low_urban_stream_pct'] * z['StreamLength'] / 1000
 
     z['CN'] = geop_result['cn']
-    z['SedPhos'] = geop_result['sed_phos']
+    z['SedPhos'] = geop_result['soilp']
     z['Area'] = [percent * area * HECTARES_PER_SQM
                  for percent in geop_result['landuse_pcts']]
 
@@ -150,7 +148,9 @@ def collect_data(geop_result, geojson):
 
     z['P'] = p_factors(z['AvSlope'], ag_lscp)
 
-    z['SedNitr'] = geop_result['soiln'] * 4.0
+    z['SedNitr'] = geop_result['soiln']
+
+    z['RecessionCoef'] = geop_result['recess_coef']
 
     return z
 
@@ -244,41 +244,38 @@ def nlcd_streams_drb(result):
 
 
 @shared_task(throws=Exception)
-def nlcd_soils(result):
+def nlcd_soil(result):
     """
     Results are expected to be in the format:
     {
-      (NLCD ID, Soil Group ID, Soil Texture ID): Count,
+      (NLCD ID, Soil Group ID): Count,
     }
 
     We calculate a number of values relying on various combinations
     of these raster datasets.
     """
     if 'error' in result:
-        raise Exception('[nlcd_soils] {}'.format(result['error']))
+        raise Exception('[nlcd_soil] {}'.format(result['error']))
 
-    ngt_count = parse(result)
+    ng_count = parse(result)
 
     # Raise exception if no NLCD values
-    if len(ngt_count.values()) == 0:
+    if len(ng_count.values()) == 0:
         raise Exception(NO_LAND_COVER)
 
     # Split combined counts into separate ones for processing
     # Reduce [(n, g, t): c] to
     n_count = {}   # [n: sum(c)]
-    ng_count = {}  # [(n, g): sum(c)]
-    nt_count = {}  # [(n, t): sum(c)]
-    for (n, g, t), count in ngt_count.iteritems():
+    ng2_count = {}  # [(n, g): sum(c)]
+    for (n, g), count in ng_count.iteritems():
         n_count[n] = count + n_count.get(n, 0)
-        nt_count[(n, t)] = count + nt_count.get((n, t), 0)
 
         # Map soil group values to usable subset
         g2 = settings.SOIL_GROUP[g]
-        ng_count[(n, g2)] = count + ng_count.get((n, g2), 0)
+        ng2_count[(n, g2)] = count + ng2_count.get((n, g2), 0)
 
     return {
-        'cn': curve_number(n_count, ng_count),
-        'sed_phos': sediment_phosphorus(nt_count),
+        'cn': curve_number(n_count, ng2_count),
         'landuse_pcts': landuse_pcts(n_count),
     }
 
@@ -319,6 +316,47 @@ def avg_awc(result):
     }
 
 
+@shared_task(throws=Exception)
+def soilp(result):
+    """
+    Get `SoilP` from MMW-Geoprocessing endpoint
+
+    Originally calculated via lookup table at Class1.vb@1.3.0:8975-8988
+    """
+    if 'error' in result:
+        raise Exception('[soilp] {}'
+                        .format(result['error']))
+
+    result = parse(result)
+
+    soilp = result.values()[0] * 1.6
+
+    return {
+        'soilp': soilp
+    }
+
+
+@shared_task(throws=Exception)
+def recess_coef(result):
+    """
+    Get `RecessCoef` from MMW-Geoprocessing endpoint
+
+    Originally a static value 0.06 Class1.vb@1.3.0:10333
+    """
+    if 'error' in result:
+        raise Exception('[recess_coef] {}'
+                        .format(result['error']))
+
+    result = parse(result)
+
+    recess_coef = result.values()[0] * -0.0015 + 0.1103
+    recess_coef = recess_coef if recess_coef >= 0 else 0.01
+
+    return {
+        'recess_coef': recess_coef
+    }
+
+
 @shared_task(throws=Exception)
 def soiln(result):
     """
@@ -332,8 +370,10 @@ def soiln(result):
 
     result = parse(result)
 
+    soiln = result.values()[0] * 4.0
+
     return {
-        'soiln': result.values()[0]
+        'soiln': soiln
     }
 
 
@@ -431,8 +471,10 @@ def nlcd_kfactor(result):
 
 def geoprocessing_chains(aoi, wkaoi, errback):
     task_defs = [
-        ('nlcd_soils',   nlcd_soils,   {'polygon': [aoi]}),
+        ('nlcd_soil',    nlcd_soil,    {'polygon': [aoi]}),
         ('soiln',        soiln,        {'polygon': [aoi]}),
+        ('soilp',        soilp,        {'polygon': [aoi]}),
+        ('recess_coef',  recess_coef,  {'polygon': [aoi]}),
         ('gwn',          gwn,          {'polygon': [aoi]}),
         ('avg_awc',      avg_awc,      {'polygon': [aoi]}),
         ('nlcd_slope',   nlcd_slope,   {'polygon': [aoi]}),

src/mmw/apps/modeling/tasks.py

@@ -69,9 +69,9 @@ def format_runoff(model_output):
 
 
 @shared_task(throws=Exception)
-def nlcd_soil_census(result):
+def nlcd_soil(result):
     if 'error' in result:
-        raise Exception('[nlcd_soil_census] {}'.format(result['error']))
+        raise Exception('[nlcd_soil] {}'.format(result['error']))
 
     dist = {}
     total_count = 0

src/mmw/apps/modeling/tests.py

@@ -128,7 +128,7 @@ def test_census(self):
                 'd:woody_wetlands': {'cell_count': 1093}
             }
         }]
-        actual = tasks.nlcd_soil_census(histogram)
+        actual = tasks.nlcd_soil(histogram)
         self.assertEqual(actual, expected)
 
 
@@ -300,7 +300,7 @@ def test_tr55_chain_doesnt_generate_censuses_if_they_exist(self):
 
         skipped_tasks = [
             'run',
-            'nlcd_soil_census'
+            'nlcd_soil'
         ]
 
         needed_tasks = [
@@ -337,7 +337,7 @@ def test_tr55_chain_doesnt_generate_aoi_census_if_it_exists_and_mods(self):
         # we still need to generate modification censuses
         needed_tasks = [
             'run',
-            'nlcd_soil_census',
+            'nlcd_soil',
             'run_tr55'
         ]
 
@@ -383,7 +383,7 @@ def test_tr55_chain_doesnt_generate_aoi_census_if_it_exists_and_no_mods(self):
 
         skipped_tasks = [
             'run',
-            'nlcd_soil_census',
+            'nlcd_soil',
         ]
 
         needed_tasks = [
@@ -443,7 +443,7 @@ def test_tr55_chain_generates_modification_censuses_if_they_are_old(self):
 
         needed_tasks = [
             'run',
-            'nlcd_soil_census',
+            'nlcd_soil',
             'run_tr55'
         ]
 
@@ -470,7 +470,7 @@ def test_tr55_chain_generates_both_censuses_if_they_are_missing(self):
 
         needed_tasks = [
             'run',
-            'nlcd_soil_census',
+            'nlcd_soil',
             'run_tr55'
         ]
 

src/mmw/apps/modeling/views.py

@@ -308,13 +308,13 @@ def _construct_tr55_job_chain(model_input, job_id):
 
         job_chain.append(tasks.run_tr55.s(censuses, aoi, model_input))
     else:
-        job_chain.append(tasks.nlcd_soil_census.s())
+        job_chain.append(tasks.nlcd_soil.s())
 
         if aoi_census and pieces:
             polygons = [m['shape']['geometry'] for m in pieces]
             geop_input = {'polygon': [json.dumps(p) for p in polygons]}
 
-            job_chain.insert(0, geoprocessing.run.s('nlcd_soil_census',
+            job_chain.insert(0, geoprocessing.run.s('nlcd_soil',
                                                     geop_input))
             job_chain.append(tasks.run_tr55.s(aoi, model_input,
                                               cached_aoi_census=aoi_census))
@@ -324,7 +324,7 @@ def _construct_tr55_job_chain(model_input, job_id):
             # Use WKAoI only if there are no pieces to modify the AoI
             wkaoi = wkaoi if not pieces else None
 
-            job_chain.insert(0, geoprocessing.run.s('nlcd_soil_census',
+            job_chain.insert(0, geoprocessing.run.s('nlcd_soil',
                                                     geop_input,
                                                     wkaoi))
             job_chain.append(tasks.run_tr55.s(aoi, model_input))

src/mmw/mmw/settings/base.py

@@ -15,7 +15,7 @@
 from layer_settings import (LAYER_GROUPS, VIZER_URLS, VIZER_IGNORE, VIZER_NAMES,
                             NHD_REGION2_PERIMETER, DRB_PERIMETER, CONUS_PERIMETER)  # NOQA
 from gwlfe_settings import (GWLFE_DEFAULTS, GWLFE_CONFIG, SOIL_GROUP, # NOQA
-                            SOILP, CURVE_NUMBER, NODATA)  # NOQA
+                            CURVE_NUMBER, NODATA)  # NOQA
 from tr55_settings import (NLCD_MAPPING, SOIL_MAPPING)
 
 # Normally you should not import ANYTHING from Django directly
@@ -431,7 +431,7 @@ def get_env_setting(setting):
                 'zoom': 0
             }
         },
-        'nlcd_soil_census': {
+        'nlcd_soil': {
             'input': {
                 'polygon': [],
                 'polygonCRS': 'LatLng',
@@ -458,20 +458,6 @@ def get_env_setting(setting):
                 'zoom': 0
             }
         },
-        'nlcd_soils': {
-            'input': {
-                'polygon': [],
-                'polygonCRS': 'LatLng',
-                'rasters': [
-                    'nlcd-2011-30m-epsg5070-512-int8',
-                    'ssurgo-hydro-groups-30m-epsg5070-512-int8',
-                    'us-ssurgo-texture-id-30m-epsg5070-512'
-                ],
-                'rasterCRS': 'ConusAlbers',
-                'operationType': 'RasterGroupedCount',
-                'zoom': 0
-            }
-        },
         'gwn': {
             'input': {
                 'polygon': [],
@@ -567,6 +553,30 @@ def get_env_setting(setting):
                 'operationType': 'RasterGroupedAverage',
                 'zoom': 0
             }
+        },
+        'soilp': {
+            'input': {
+                'polygon': [],
+                'polygonCRS': 'LatLng',
+                'rasters': [],
+                'targetRaster': 'soilpallland2-epsg5070',
+                'pixelIsArea': True,
+                'rasterCRS': 'ConusAlbers',
+                'operationType': 'RasterGroupedAverage',
+                'zoom': 0
+            }
+        },
+        'recess_coef': {
+            'input': {
+                'polygon': [],
+                'polygonCRS': 'LatLng',
+                'rasters': [],
+                'targetRaster': 'bfi48grd-epsg5070',
+                'pixelIsArea': True,
+                'rasterCRS': 'ConusAlbers',
+                'operationType': 'RasterGroupedAverage',
+                'zoom': 0
+            }
         }
     }
 }