Sed_export calculation with drainage raster instead of stream raster (InVEST SDR)

[Marvin]

Hi everyone,

i have a very detailed raster of streams and lakes, very accurate compared to the stream.tif that is calculated by the SDR model based on the relief.
I would like to calculate the sed_export.tif raster based on this accurate raster but exclude the stream.tif that is usually used in the InVEST SDR model. An optional drainage raster can be add to the model and is combined with the stream.tif raster to stream_and_drainage.tif. So i can use this option to add my own raster of streams and lakes. Is it possible to change source code so only the drainage raster determines where soil is exported to a stream?

Maybe changes have to be done here:

def add_drainage_op(stream, drainage):
“”“Add drainage mask to stream layer.”“”
return numpy.where(drainage == 1, 1, stream)

stream_nodata = pygeoprocessing.get_raster_info(stream_path)['nodata'][0]
pygeoprocessing.raster_calculator(
    [(path, 1) for path in [stream_path, drainage_path]], add_drainage_op,
    out_stream_and_drainage_path, gdal.GDT_Byte, stream_nodata)

I would be pleased to receive suggestions.

Best,
Marvin

[esoth]

Hi @Marvin,
Yes, you could definitely modify the source code to do that! You’re right that add_drainage_op would be the place to modify. The statement numpy.where(drainage == 1, 1, stream) is equivalent to drainage | stream. It’s doing a logical OR on the two rasters so that the output has a 1 where either drainage or stream has a 1.

The quickest way to make it only consider your drainage raster would be to replace the line return numpy.where(drainage == 1, 1, stream) with return drainage.
This is kind of redundant; really you could remove the whole _add_drainage task, but then you’d have to modify the execute function which is more complicated.
This sounds like a useful feature, though. I’ll see if it would be possible to modify the UI, so in a future release of InVEST, you could choose to use either the known drainage raster, the calculated stream raster, or both.

I hope that helps. Let us know if you have further questions about modifying and running the source code!

[Marvin]

Hi @esoth ,
thanks a lot for the quick reply! I completely agree, this would be a nice feature to have. Especially on arable land natural streams are regulated and partly in the subsurface. In this cases the natural streams calculated based on the relief are not sufficient for detailed calculations.
I already ran the model with the modified source code using > return drainage.

As a test I used a drainage grid that only has pixel values = 0. So there are no rivers within the analysed area. Nevertheless, I get a sed_export raster with pixel values > 0. In the guidelines sed_export is defined as “The total amount of sediment exported from each pixel that reaches the stream.” Why is it possible that an export takes place even though there are no rivers?

Thanks for your help!

Marvin

[esoth]

Hi @Marvin,
The sediment export algorithm depends on the sediment delivery ratio (SDR), which depends on the path and distance to the nearest stream. When there are no streams, it doesn’t make sense to calculate the SDR in this way. Really, the SDR and sediment export rasters should be all nodata when there aren’t any streams.

I’m not sure exactly where the sed_export values are coming from in this case, but the distance values aren’t “real”, so the SDR and sed_export values aren’t “real” either. If you have at least one stream pixel, the values should make more sense.

[Marvin]

Hi @esoth,

i tried a drainage raster that includes one “stripe” with pixel values = 1. As a result i get a sed_export.tif which has higher erosion values near the drainage stripe compared to calculations were i used a drainage raster with pixel values = 0. But still there is sediment export in the results far away from the stripe. I overlapped the drainage raster with the stripe (darker blue), the stream.tif (lighter blue) and the sed_export.tif. For me it seems that the stream.tif still has an influence on the sed_export.tif (see figure). Is there an other part of source code that has to be changed?

Thanks,
Marvin

[esoth]

Hi @Marvin,
Are you running python setup.py install after making changes to the source code? That step is necessary to update the natcap.invest package with the changes you’ve made.

[Marvin]

Hi @esoth ,
my approach (possibly a little cumbersome):
I have created a python environment in Anaconda.
I use an unpacked whl file of InVEST 3.9 where I change the source code. Then I pack the file with wheel pack. Then I uninstall the old version with pip uninstall natcap.invest
and install the new version with pip install. Is this a possible solution?

I guess there is an easier way?

Thanks,
Marvin

[esoth]

Hi @Marvin,
That sounds like it should work, but I’m not certain. The workflow I’d recommend is:

1. Fork the source code from natcap/invest
2. Make your changes and commit them
3. Run python setup.py install
4. Run pip list or conda list and double check that the natcap.invest version number is the same as the version number that was printed out at the end of python setup.py install

Some benefits of doing it this way:

• You’re creating the package in the same way that we do for versions we release. While I don’t think it matters for this change you’re making now, the setup.py script can make important changes/settings that don’t get called if you only unpack/repack the wheel.
• If you ever want to modify one of the Cython modules, python setup.py install is necessary to compile it. Otherwise you’d have to directly modify the .so file in the wheel (not recommended).
• You get the benefits of version control and can compare directly to the original 3.9.0 version and our latest changes.

[Marvin]

Hi @esoth ,

i downloaded the current code from the link you sent me. When i try python setup.py install i get the following error:

(INVEST38) C:\invest-main>python setup.py install
Traceback (most recent call last):
File “setup.py”, line 40, in
setup(
File “C:\Users\melzerm\Anaconda3\envs\INVEST38\lib\site-packages\setuptools_init_.py”, line 153, in setup
return distutils.core.setup(**attrs)
File “C:\Users\melzerm\Anaconda3\envs\INVEST38\lib\distutils\core.py”, line 108, in setup
_setup_distribution = dist = klass(attrs)
File “C:\Users\melzerm\Anaconda3\envs\INVEST38\lib\site-packages\setuptools\dist.py”, line 434, in init
_Distribution.init(self, {
File “C:\Users\melzerm\Anaconda3\envs\INVEST38\lib\distutils\dist.py”, line 292, in init
self.finalize_options()
File “C:\Users\melzerm\Anaconda3\envs\INVEST38\lib\site-packages\setuptools\dist.py”, line 743, in finalize_options
ep(self)
File “C:\Users\melzerm\Anaconda3\envs\INVEST38\lib\site-packages\setuptools\dist.py”, line 750, in _finalize_setup_keywords
ep.load()(self, ep.name, value)
File “c:\invest-main.eggs\setuptools_scm-6.0.0-py3.8.egg\setuptools_scm\integration.py”, line 24, in version_keyword
dist.metadata.version = get_version(config)
File "c:\invest-main.eggs\setuptools_scm-6.0.0-py3.8.egg\setuptools_scm_init.py", line 173, in _get_version
parsed_version = do_parse(config)
File "c:\invest-main.eggs\setuptools_scm-6.0.0-py3.8.egg\setuptools_scm_init.py", line 134, in _do_parse
raise LookupError(
LookupError: setuptools-scm was unable to detect version for ‘C:\invest-main’.

Make sure you’re either building from a fully intact git repository or PyPI tarballs. Most other sources (such as GitHub’s tarballs, a git checkout without the .git folder) don’t contain the necessary metadata and will not work.

For example, if you’re using pip, instead of https://github.com/user/proj/archive/master.zip use git+https://github.com/user/proj.git#egg=proj

Do you know the reason for this error? Regardless of this approach i think my solution works fine for this small source code modification. But still i am searching for reasons why sed_export contains pixel values >0 in areas without drainage or streams.

Thanks,
Marvin

[dcdenu4]

Hey @Marvin ,

Sorry for the delay, we only have limited support for source modifications. If I’m following the discussion so far you’re wondering why you are still seeing sed_export values where there are no streams.

This is because the SDR model is calculating a flow accumulation layer that it uses for the ls_factor which is then used in the rkls calculation. flow accumulation is dependent on the DEM. The stream layer is derived from the flow accumulation layer and is largely used in the model to mask out values if they happen to be “in” a stream. So why you could certainly bump the Threshold flow accumulation value so that you would have a stream layer of all 0s, you would still end up with sed_export values. However, by setting the threshold flow accumulation value to extremely high, you could then add your own drainage layer with a more accurate stream representation. However, this layer only indicates which pixels drain to the watershed.

If you really wanted to alter the model to only look at export using your own drainage layer, you would need to create your own flow accumulation raster and substitute that into the source as well. Sorry this is not more straight forward.

Cheers,

Doug

[Marvin]

Hey @dcdenu4 ,

thanks a lot for the explanation! So how should this own flow accumulation raster look like? Can it be the same raster as the one i use for the drainage layer? Have you any suggestions how to substitute the flow accumulation raster to the model? Do i have to make further changes on the source code or just copy the raster into the output folder were the regular flow accumulation raster is saved?

Thanks,

Marvin

[dcdenu4]

Hi @Marvin,

Flow Accumulation is a raster where: for current pixel p_i the value of p_i is the number of upstream pixels that enter into p_i. This helps us compute the contributing area by multiplying by the cell_area. This is used in the LS equation as shown in the Users Guide here: https://storage.googleapis.com/releases.naturalcapitalproject.org/invest-userguide/latest/sdr.html#annual-soil-loss

Can it be the same raster as the one i use for the drainage layer?

No, these are separate and have different meanings. The drainage layer is just an etching of which pixels exit downstream and out of the watershed. And also, which pixels should not have SDR related outputs because it’s a waterway.

Have you any suggestions how to substitute the flow accumulation raster to the model?

I would first suggest using your drainage layer as input and setting the threshold value to unreasonably high so that you get a stream layer of all 0s. I would take a look at these results and see if how they fit your data, and to get a sense of what the model and what the drainage layer is doing. It might also be worth playing around with the flow accumulation threshold to see how the stream orders can change.

From there you could do a few things off the top of my head:

• Burn your drainage layer into the DEM by slightly raising DEM values around your drainage layer. Then use that DEM as input into the model and see if the stream layer then aligns better with your expected drainage layer. Note, there are consequences to this because slope is a big factor in the model and by altering the DEM you are also altering slope values.
• Take the SDR produced Flow Accumulation and manually set values to 0 where you know pixels should not be contributing. Then edit source to use that raster specifically.
• Consider that this SDR model does not suit your needs. This is a general purpose model that is based on the USLE equation. Read the Limitation section of the User’s Guide: https://storage.googleapis.com/releases.naturalcapitalproject.org/invest-userguide/latest/sdr.html#limitations

Do i have to make further changes on the source code or just copy the raster into the output folder were the regular flow accumulation raster is saved?

It’s complicated, but I think it would be easiest to just use the path of the modified flow accumulation raster in the code directly, commenting out the code that calls flow accumulation, and replacing the path variable that the model expects with the new modified path.

This is confusing and complicated because it is stretching the limits on what this model is intended to do and the problem it was designed to solve.

Doug

[Marvin]

Thanks a lot @dcdenu4 for this detailed answer! I will try some of your suggestions. Setting the threshold value high seems to be a satisfactory solution for my aims.

Thanks,
Marvin