# Sed\_export calculation with drainage raster instead of stream raster (InVEST SDR)

**URL:** <https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757>\
**Category:** Software Support\
**Tags:** sdr\
**Created:** [March 10, 2021, 1:35pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757 "2021-03-10T13:35:56Z")\
**Posts on this page:** 13\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marvin](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/marvin/32/1214_2.png) [@Marvin](https://community.naturalcapitalalliance.org/u/Marvin)\
**Post date:** [March 10, 2021, 1:35pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/1 "2021-03-10T13:35:56Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![esoth](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@esoth](https://community.naturalcapitalalliance.org/u/esoth)\
**Post date:** [March 10, 2021, 6:45pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/2 "2021-03-10T18:45:33Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![Marvin](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/marvin/32/1214_2.png) [@Marvin](https://community.naturalcapitalalliance.org/u/Marvin)\
**Post date:** [March 11, 2021, 7:28am UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/3 "2021-03-11T07:28:21Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![esoth](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@esoth](https://community.naturalcapitalalliance.org/u/esoth)\
**Post date:** [March 11, 2021, 7:14pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/4 "2021-03-11T19:14:49Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Marvin](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/marvin/32/1214_2.png) [@Marvin](https://community.naturalcapitalalliance.org/u/Marvin)\
**Post date:** [March 12, 2021, 8:02am UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/5 "2021-03-12T08:02:18Z")

</div>

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?

 ![Drainage_and_stream](https://canada1.discourse-cdn.com/flex031/uploads/naturalcapitalproject/original/2X/d/d54f3a2ab21e4568315471a9239d3997eeb6f7bc.jpeg)

Thanks,  
Marvin

---

<div class="post-metadata">

**Author:** ![esoth](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@esoth](https://community.naturalcapitalalliance.org/u/esoth)\
**Post date:** [March 12, 2021, 5:09pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/6 "2021-03-12T17:09:55Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Marvin](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/marvin/32/1214_2.png) [@Marvin](https://community.naturalcapitalalliance.org/u/Marvin)\
**Post date:** [March 15, 2021, 8:04am UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/7 "2021-03-15T08:04:14Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![esoth](https://avatars.discourse-cdn.com/v4/letter/e/9d8465/32.png) [@esoth](https://community.naturalcapitalalliance.org/u/esoth)\
**Post date:** [March 15, 2021, 4:31pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/8 "2021-03-15T16:31:34Z")

</div>

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](https://github.com/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.

---

<div class="post-metadata">

**Author:** ![Marvin](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/marvin/32/1214_2.png) [@Marvin](https://community.naturalcapitalalliance.org/u/Marvin)\
**Post date:** [March 17, 2021, 8:45am UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/9 "2021-03-17T08:45:38Z")

</div>

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](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

---

<div class="post-metadata">

**Author:** ![dcdenu4](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/dcdenu4/32/544_2.png) [@dcdenu4](https://community.naturalcapitalalliance.org/u/dcdenu4)\
**Post date:** [March 23, 2021, 1:34pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/10 "2021-03-23T13:34:51Z")

</div>

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 `0`s, 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

---

<div class="post-metadata">

**Author:** ![Marvin](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/marvin/32/1214_2.png) [@Marvin](https://community.naturalcapitalalliance.org/u/Marvin)\
**Post date:** [March 23, 2021, 2:25pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/11 "2021-03-23T14:25:25Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![dcdenu4](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/dcdenu4/32/544_2.png) [@dcdenu4](https://community.naturalcapitalalliance.org/u/dcdenu4)\
**Post date:** [March 23, 2021, 5:33pm UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/12 "2021-03-23T17:33:45Z")

</div>

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](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](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

---

<div class="post-metadata">

**Author:** ![Marvin](https://yyz1.discourse-cdn.com/flex031/user_avatar/community.naturalcapitalalliance.org/marvin/32/1214_2.png) [@Marvin](https://community.naturalcapitalalliance.org/u/Marvin)\
**Post date:** [April 30, 2021, 6:11am UTC](https://community.naturalcapitalalliance.org/t/sed-export-calculation-with-drainage-raster-instead-of-stream-raster-invest-sdr/1757/13 "2021-04-30T06:11:10Z")

</div>

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
