re #BITBUCKET-3626 removing regression function because I don't know what it'l be yet and removing the corresponding test in suite.
This commit is contained in:
Rich Sharp
parent
b538e2657c
commit
1f70f5d89b
|
|
@ -1,221 +0,0 @@
|
|||
"""InVEST Crop Production Regression Model."""
|
||||
import collections
|
||||
import re
|
||||
import os
|
||||
import logging
|
||||
|
||||
import numpy
|
||||
from osgeo import gdal
|
||||
from osgeo import osr
|
||||
from osgeo import ogr
|
||||
import pygeoprocessing
|
||||
|
||||
from . import utils
|
||||
|
||||
logging.basicConfig(format='%(asctime)s %(name)-20s %(levelname)-8s \
|
||||
%(message)s', level=logging.DEBUG, datefmt='%m/%d/%Y %H:%M:%S ')
|
||||
|
||||
LOGGER = logging.getLogger('natcap.invest.crop_production_regression')
|
||||
|
||||
_INTERMEDIATE_OUTPUT_DIR = 'intermediate_output'
|
||||
|
||||
_YIELD_PERCENTILE_FIELD_PATTERN = 'yield_([^_]+)'
|
||||
_GLOBAL_OBSERVED_YIELD_FILE_PATTERN = os.path.join(
|
||||
'observed_yield', '%s_yield_map.tif') # crop_name
|
||||
_EXTENDED_CLIMATE_BIN_FILE_PATTERN = os.path.join(
|
||||
'extended_climate_bin_maps', 'extendedclimatebins%s.tif') # crop_name
|
||||
_CLIMATE_PERCENTILE_TABLE_PATTERN = os.path.join(
|
||||
'climate_percentile_yield_tables',
|
||||
'%s_percentile_yield_table.csv') # crop_name
|
||||
|
||||
# crop_name, yield_percentile_id
|
||||
_INTERPOLATED_YIELD_PERCENTILE_FILE_PATTERN = os.path.join(
|
||||
_INTERMEDIATE_OUTPUT_DIR, '%s_%s_interpolated_yield%s.tif')
|
||||
|
||||
# crop_name, file_suffix
|
||||
_CLIPPED_CLIMATE_BIN_FILE_PATTERN = os.path.join(
|
||||
_INTERMEDIATE_OUTPUT_DIR,
|
||||
'clipped_%s_climate_bin_map%s.tif')
|
||||
|
||||
# crop_name, yield_percentile_id, file_suffix
|
||||
_COARSE_YIELD_PERCENTILE_FILE_PATTERN = os.path.join(
|
||||
_INTERMEDIATE_OUTPUT_DIR, '%s_%s_coarse_yield%s.tif')
|
||||
|
||||
# crop_name, yield_percentile_id, file_suffix
|
||||
_PERCENTILE_CROP_PRODUCTION_FILE_PATTERN = os.path.join(
|
||||
'.', '%s_%s_production%s.tif')
|
||||
|
||||
# crop_name, file_suffix
|
||||
_CLIPPED_OBSERVED_YIELD_FILE_PATTERN = os.path.join(
|
||||
_INTERMEDIATE_OUTPUT_DIR, '%s_clipped_observed_yield%s.tif')
|
||||
|
||||
# crop_name, file_suffix
|
||||
_ZEROED_OBSERVED_YIELD_FILE_PATTERN = os.path.join(
|
||||
_INTERMEDIATE_OUTPUT_DIR, '%s_zeroed_observed_yield%s.tif')
|
||||
|
||||
# crop_name, file_suffix
|
||||
_INTERPOLATED_OBSERVED_YIELD_FILE_PATTERN = os.path.join(
|
||||
_INTERMEDIATE_OUTPUT_DIR, '%s_interpolated_observed_yield%s.tif')
|
||||
|
||||
# crop_name, file_suffix
|
||||
_OBSERVED_PRODUCTION_FILE_PATTERN = os.path.join(
|
||||
'.', '%s_observed_production%s.tif')
|
||||
|
||||
# file_suffix
|
||||
_AGGREGATE_VECTOR_FILE_PATTERN = os.path.join(
|
||||
_INTERMEDIATE_OUTPUT_DIR, 'aggrgate_vector%s.shp')
|
||||
|
||||
# file_suffix
|
||||
_AGGREGATE_TABLE_FILE_PATTERN = os.path.join(
|
||||
'.', 'aggregate_results%s.csv')
|
||||
|
||||
_EXPECTED_NUTRIENT_TABLE_HEADERS = [
|
||||
'Protein', 'Lipid', 'Energy', 'Ca', 'Fe', 'Mg', 'Ph', 'K', 'Na', 'Zn',
|
||||
'Cu', 'Fl', 'Mn', 'Se', 'VitA', 'betaC', 'alphaC', 'VitE', 'Crypto',
|
||||
'Lycopene', 'Lutein', 'betaT', 'gammaT', 'deltaT', 'VitC', 'Thiamin',
|
||||
'Riboflavin', 'Niacin', 'Pantothenic', 'VitB6', 'Folate', 'VitB12',
|
||||
'VitK']
|
||||
_EXPECTED_LUCODE_TABLE_HEADER = 'lucode'
|
||||
_NODATA_CLIMATE_BIN = 255
|
||||
_NODATA_YIELD = -1.0
|
||||
|
||||
|
||||
def execute(args):
|
||||
"""Crop Production Regression Model.
|
||||
|
||||
This model will take a landcover, information about fertilization of
|
||||
K, N, and Pot, and irrigation information to produce modeled ag
|
||||
production.
|
||||
|
||||
[say equations here and cite reference?]
|
||||
|
||||
Parameters:
|
||||
args['workspace_dir'] (string): output directory for intermediate,
|
||||
temporary, and final files
|
||||
args['results_suffix'] (string): (optional) string to append to any
|
||||
output file names
|
||||
args['landcover_raster_path'] (string): path to landcover raster
|
||||
args['landcover_to_crop_table_path'] (string): path to a table that
|
||||
converts landcover types to crop names that has two headers:
|
||||
* lucode: integer value corresponding to a landcover code in
|
||||
`args['landcover_raster_path']`.
|
||||
* crop_name: a string that must match one of the crops in
|
||||
args['model_data_path']/climate_bin_maps/[cropname]_*
|
||||
A ValueError is raised if strings don't match.
|
||||
args['k_raster_path'] (string): path to K fertilization rates (kg/Ha).
|
||||
args['n_raster_path'] (string): path to N fertilization rates (kg/Ha).
|
||||
args['pot_raster_path'] (string): path to Pot fertilization rates (
|
||||
kg/Ha).
|
||||
args['irrigation_raster_path'] (string): path to irrigation mask,
|
||||
pixel value of 1 indicates irrigation present while 0 does not.
|
||||
args['aggregate_polygon_path'] (string): path to polygon shapefile
|
||||
that will be used to aggregate crop yields and total nutrient
|
||||
value. (optional, if value is None, then skipped)
|
||||
args['aggregate_polygon_id'] (string): if
|
||||
args['aggregate_polygon_path'] is provided, then this value is a
|
||||
id field in that vector that will be used to index the final
|
||||
aggregate results.
|
||||
args['model_data_path'] (string): path to the InVEST Crop Production
|
||||
global data directory. This model expects that the following
|
||||
directories are subdirectories of this path
|
||||
* climate_bin_maps (contains [cropname]_climate_bin.tif files)
|
||||
* climate_percentile_yield (contains
|
||||
[cropname]_percentile_yield_table.csv files)
|
||||
|
||||
Returns:
|
||||
None.
|
||||
"""
|
||||
file_suffix = utils.make_suffix_string(args, 'results_suffix')
|
||||
|
||||
output_dir = os.path.join(args['workspace_dir'])
|
||||
utils.make_directories([
|
||||
output_dir, os.path.join(output_dir, _INTERMEDIATE_OUTPUT_DIR)])
|
||||
|
||||
landcover_raster_info = pygeoprocessing.get_raster_info(
|
||||
args['landcover_raster_path'])
|
||||
pixel_area_ha = numpy.product([
|
||||
abs(x) for x in landcover_raster_info['pixel_size']]) / 10000.0
|
||||
landcover_nodata = landcover_raster_info['nodata'][0]
|
||||
|
||||
crop_to_landcover_table = utils.build_lookup_from_csv(
|
||||
args['landcover_to_crop_table_path'], 'crop_name', to_lower=True,
|
||||
numerical_cast=True)
|
||||
|
||||
crop_lucodes = [
|
||||
x[_EXPECTED_LUCODE_TABLE_HEADER]
|
||||
for x in crop_to_landcover_table.itervalues()]
|
||||
|
||||
LOGGER.info(
|
||||
"Calculating total land area and warning if the landcover raster "
|
||||
"is missing lucodes")
|
||||
unique_lucodes = numpy.array([])
|
||||
total_area = 0.0
|
||||
for _, lu_band_data in pygeoprocessing.iterblocks(
|
||||
args['landcover_raster_path']):
|
||||
unique_block = numpy.unique(lu_band_data)
|
||||
unique_lucodes = numpy.unique(numpy.concatenate(
|
||||
(unique_lucodes, unique_block)))
|
||||
total_area += numpy.count_nonzero((lu_band_data != _NODATA_YIELD))
|
||||
|
||||
missing_lucodes = set(crop_lucodes).difference(
|
||||
set(unique_lucodes))
|
||||
if len(missing_lucodes) > 0:
|
||||
LOGGER.warn(
|
||||
"The following lucodes are in the landcover to crop table but "
|
||||
"aren't in the landcover raster: %s", missing_lucodes)
|
||||
|
||||
# Calculate lat/lng bounding box for landcover map
|
||||
wgs84srs = osr.SpatialReference()
|
||||
wgs84srs.ImportFromEPSG(4326) # EPSG4326 is WGS84 lat/lng
|
||||
landcover_wgs84_bounding_box = pygeoprocessing.transform_bounding_box(
|
||||
landcover_raster_info['bounding_box'],
|
||||
landcover_raster_info['projection'], wgs84srs.ExportToWkt(),
|
||||
edge_samples=11)
|
||||
|
||||
crop_lucode = None
|
||||
observed_yield_nodata = None
|
||||
crop_area = collections.defaultdict(float)
|
||||
for crop_name in crop_to_landcover_table:
|
||||
crop_lucode = crop_to_landcover_table[crop_name][
|
||||
_EXPECTED_LUCODE_TABLE_HEADER]
|
||||
print crop_name, crop_lucode
|
||||
crop_climate_bin_raster_path = os.path.join(
|
||||
args['model_data_path'],
|
||||
_EXTENDED_CLIMATE_BIN_FILE_PATTERN % crop_name)
|
||||
if not os.path.exists(crop_climate_bin_raster_path):
|
||||
raise ValueError(
|
||||
"Expected climate bin map called %s for crop %s "
|
||||
"specified in %s", crop_climate_bin_raster_path, crop_name,
|
||||
args['landcover_to_crop_table_path'])
|
||||
if not os.path.exists(crop_climate_bin_raster_path):
|
||||
raise ValueError(
|
||||
"Expected climate bin map called %s for crop %s "
|
||||
"specified in %s", crop_climate_bin_raster_path, crop_name,
|
||||
args['landcover_to_crop_table_path'])
|
||||
|
||||
clipped_climate_bin_raster_path = os.path.join(
|
||||
output_dir, _CLIPPED_CLIMATE_BIN_FILE_PATTERN % (
|
||||
crop_name, file_suffix))
|
||||
|
||||
crop_climate_bin_raster_info = pygeoprocessing.get_raster_info(
|
||||
crop_climate_bin_raster_path)
|
||||
pygeoprocessing.warp_raster(
|
||||
crop_climate_bin_raster_path,
|
||||
crop_climate_bin_raster_info['pixel_size'],
|
||||
clipped_climate_bin_raster_path, 'nearest',
|
||||
target_bb=landcover_wgs84_bounding_box)
|
||||
|
||||
climate_percentile_yield_table_path = os.path.join(
|
||||
args['model_data_path'],
|
||||
_CLIMATE_PERCENTILE_TABLE_PATTERN % crop_name)
|
||||
crop_climate_percentile_table = utils.build_lookup_from_csv(
|
||||
climate_percentile_yield_table_path, 'climate_bin',
|
||||
to_lower=True, numerical_cast=True)
|
||||
|
||||
yield_percentile_headers = [
|
||||
x for x in crop_climate_percentile_table.itervalues().next()
|
||||
if x != 'climate_bin']
|
||||
|
||||
clipped_climate_bin_raster_path_info = (
|
||||
pygeoprocessing.get_raster_info(
|
||||
clipped_climate_bin_raster_path))
|
||||
|
|
@ -1,10 +1,10 @@
|
|||
"""Module for Regression Testing the InVEST Crop Production models."""
|
||||
import numpy
|
||||
import unittest
|
||||
import tempfile
|
||||
import shutil
|
||||
import os
|
||||
|
||||
import numpy
|
||||
import pygeoprocessing.testing
|
||||
from pygeoprocessing.testing import scm
|
||||
|
||||
|
|
@ -106,30 +106,3 @@ class CropProductionTests(unittest.TestCase):
|
|||
|
||||
with self.assertRaises(ValueError):
|
||||
crop_production_percentile.execute(args)
|
||||
|
||||
|
||||
@unittest.skip("Still in development.")
|
||||
@scm.skip_if_data_missing(SAMPLE_DATA_PATH)
|
||||
@scm.skip_if_data_missing(MODEL_DATA_PATH)
|
||||
def test_crop_production_regression(self):
|
||||
"""Crop Production: test crop production."""
|
||||
from natcap.invest import crop_production_regression
|
||||
|
||||
args = {
|
||||
'workspace_dir': self.workspace_dir,
|
||||
'results_suffix': '',
|
||||
'landcover_raster_path': os.path.join(
|
||||
SAMPLE_DATA_PATH, 'landcover.tif'),
|
||||
'landcover_to_crop_table_path': os.path.join(
|
||||
SAMPLE_DATA_PATH, 'landcover_to_crop_table.csv'),
|
||||
'k_raster_path': os.path.join(
|
||||
SAMPLE_DATA_PATH, 'k.tif'),
|
||||
'n_raster_path': os.path.join(
|
||||
SAMPLE_DATA_PATH, 'n.tif'),
|
||||
'pot_raster_path': os.path.join(
|
||||
SAMPLE_DATA_PATH, 'pot.tif'),
|
||||
'irrigation_raster_path': os.path.join(
|
||||
SAMPLE_DATA_PATH, 'irrigation.tif'),
|
||||
'model_data_path': MODEL_DATA_PATH
|
||||
}
|
||||
crop_production_regression.execute(args)
|
||||
|
|
|
|||
Loading…
Reference in New Issue