Input datasets for SDR (data sources & pre-processing)

[ndmetherall]

Dear @swolny - many thanks again for the technical support and software troubleshooting advice.

I have taken your advice and followed up on the following:

1. I have tried to use a different filled sinks version of the DEM
2. I have also stopped using my work laptop which for some reason automatically sinks to a dropbox
3. I have got a new license to install on my desktop computer which is not automatically saving to a dropbox
4. I am now using ArcMAP rather than ArcGIS Pro - this is allowing me to save in .TIFF format

Shapefile with revised attribute table Uploading: 5_W2.zip…

7. Finally my LULC raster and my erodibility raster don’t fully cover the extent of the Southernmost island in this screenshot below. Can the model still run for the area that fits in the extent or will this stop the whole model from running?

---

The next day I tried to run the model again and this time it ran for some reason. According to the logfile attached there are no errors this time.
InVEST-Sediment-Delivery-Ratio-Model-(SDR)-log-2021-07-24–10_28_21.txt (1.1 MB)

However, the outputs don’t look quite right either see below:

Watersheds vector seems to be ok:

However the raster outputs do not seem right judging from the values ranges - sediment deposition is negative :

USLE results don’t provide much insight as shown here. Although the southernmost island is entirely blacked out.

A zoomed in version shows the USLE maybe showing deposition on the coastal ocean?

Sediment deposition just gives a high value for all ocean and a low value for all land
Processing: image.png…

All advice is appreciated on how to troubleshoot this is appreciated. Looking forward to discussing.

Thanks again

[swolny]

@ndmetherall

7. Finally my LULC raster and my erodibility raster don’t fully cover the extent of the Southernmost island in this screenshot below. Can the model still run for the area that fits in the extent or will this stop the whole model from running?

If any of your rasters have NoData values for some set of pixels, typically the model will run, but the result will be NoData in those pixels. We do of course recommend having all of the inputs cover the whole study area, otherwise your output will be incomplete in those watersheds with missing data.

@ndmetherall

However the raster outputs do not seem right judging from the values ranges - sediment deposition is negative :

Sediment deposition should definitely not have those huge negative values. Where do those pixels occur? Do you see a pattern with the areas of missing data, shoreline or anything else?

@ndmetherall

Sediment deposition just gives a high value for all ocean and a low value for all land

I’m confused why the model is producing results on the ocean at all (which is not obvious from the map; not knowing this area it looks like a bunch of landmasses to me) - do all of your inputs cover the ocean too? (I’d be surprised if soil and land use data do). Since these are islands, and presumably a lot of the watersheds flow to the sea, the data really should cut off at the shore, otherwise the model will keep applying its land-based calculations to the marine environment, which is not appropriate.

@ndmetherall

A zoomed in version shows the USLE maybe showing deposition on the coastal ocean?

It’s hard to parse this map - are the watersheds covering the USLE raster? Or are USLE results only being produced for those black coastal pixels? USLE does not indicate deposition (sed_deposition does that), it indicates the source of erosion, before it moves down the landscape.

~ Stacie

[ndmetherall]

1. Thanks for this info about the extent. I will expand the extent once I can complete the model run.

2. The sediment deposition seems to correlate mostly to the land use and DEM layer inputs. The Solomon Islands are mostly dense forest with some urban expansion. The patterns of missing data are on the ocean. The trends from ridge to shoreline just shows a uniform low value for sediment deposition.

Sediment deposition

raster legend

Overlaying DEM and LULC over sediment deposition

3. My mistake. There is no value for the ocean for sed deposition.

4. Sediment retention, sediment export and rkls all show similar results of deposition on the coastline.
   

   
   - however, I am unsure if this is accurate

5. These are our USLE results. They seem to only be produced for the coastal pixels as you have described. So does this mean that soil is not being lost from the forested landuse class upstream?
   

   

I will try to revisit the input datasets and see if I can share them to help diagnose / parse the map.

Thank you for the support.

[ndmetherall]

Trying to share the input files in zip form:

2. clipped erositivity 2_R.tif.zip (137.1 KB)
3. clipped erodibility 3_K.tif.zip (301.6 KB)
4. clipped LULC 4_LULC.tif.vat.zip (16.4 KB)
5. clipped biophysical table SI_joined_table.csv (339 Bytes)

[ndmetherall]

I have had to upload these to Google Drive in order to share due to the file sizes. Please let me know if there are any issues accessing.

1. DEM
   1_DEM.zip - Google Drive

2. Watersheds generated through DelineatIt
   5_W2.zip - Google Drive

Looking forward to discussing.

Kind regards.

[swolny]

Does the blue that’s covering most of the island have the large negative value, or are the values more reasonable? When the value ramp has a huge range like that, it will certainly make symbology difficult, and skewed.

It’s not making sense to me that the sediment deposition output covers the whole island, but USLE only covers the coast. All output layers should cover the whole area that you’ve provided inputs for. Something is very strange there.

I’ll check out and try running SDR with your inputs, thanks for uploading them. But, can you re-try uploading the watersheds? I can’t get 5_W2.shp to open in either Arc or Q, they both say it’s an unrecognized format. Thanks.

~ Stacie

[ndmetherall]

Hi @swolny

1. The blue with the large negative values cover the entire land mass. The values do not seem to be very reasonable - just a flat single value aligning with most of the land dataset.

2. I am trying to upload the watersheds again. I seem to be having the same strange issue when trying to redownload the 5_W2.shp file. It is strange for this to happen as it seems to work on the computer but maybe has been altered when uploaded to google docs somehow.

3. I hope to reupload today. I generated the 5_W2.shp using the delineatIt module which took a few hours to run - using the same DEM you have there as the input to the model. So I am hoping I will not have to rerun this but will let you know.

Best wishes.

[ndmetherall]

Dear @swolny

Please find watersheds attached from a previous DelineateIt run - the most recent run was done on a different computer.

1. shape file alone
   SI_watersheds shp.zip (204.2 KB)

2. shape file and other files: (upload://wnIzeE8gkwWsuDhLxFocI9di4OD.zip) (204.2 KB)
   SI_watersheds pkg.zip (209.8 KB)

Thanks for the help

[swolny]

Hey @ndmetherall, I’m sorry to report that I still can’t use that shapefile, neither ArcGIS nor InVEST will accept it. No idea what the problem is, but there is a strange “.fix” file, which it’s probably ignoring, but it’s also missing an .shx file, which is required. For the helluvit I tried changing .fix to .shx, but that didn’t help. How are you packaging up the shapefile? If you haven’t already, try making a new folder, and use ArcCatalog to copy the working shapefile into that folder (it will automatically include all necessary files), then zip it up. If that particular method is producing this problem, that’s very surprising.

~ Stacie

[ndmetherall]

I appreciate your patience @swolny - apologies for the errors on my side with packaging this file. I am going to try another zip with .shx here:

I have two versions:

1. one where I have just included the .shx link and all files in a zap “5_WS”
2. one where I have followed your instructions exactly, making a new folder and using ArcCatalog to copy the working shapefile into that folder. Just checking because when I do this, the options in ArcCatalog are ‘export’ or ‘copy’ and ‘export’ seems to be the one which allows me to save another copy of the feature class on my PC. The result of this export is attached below.

You can find these two zip files in the google drive folder below:
https://drive.google.com/drive/folders/1KuszKmLk48DfMPwvz87lH3NV9ZecE_O3?usp=sharing

Please let me know if this approach works.

Best wishes.

[swolny]

Thanks, those latest shapefiles both worked. So I found two problems.

One, your K/erodibility raster has NoData over most of the landmasses - this is probably why your RKLS and USLE outputs only appear along the coast, and of course this will trickle down to subsequent calculations.

Second, I ran the model anyway, and saw that no streams are being created with the threshold flow accumulation value of 1000 that was listed as being used in one of your log files - look at stream.tif. SDR moves erosion downslope until it reaches a stream, and if no streams are created, it can’t calculate export. So you’ll need to try different threshold flow accumulation values until streams are created that look reasonably close to reality.

Try fixing those two things and see what happens.

~ Stacie

[ndmetherall]

Dear @swolny - many thanks indeed for diagnosing the issues. I can work with this.

1. So firstly in terms of the erodibility raster, I will just try to redo this layer from scratch from soil grids. There may have been an issue with working on a computer automatically linked to a dropbox which you highlighted as a potential issue (thank you). Can you think of any other reasons why this might be occurring and what I might need to do differently to generate the erodibility input?

2. I will try to rerun the SDR model with multiple different flow accumulation values (i.e. 500, 250, 100 etc…). I am guessing that lower numbers will be more likely to generate streams?

I will try again and see what happens.

Many thanks again for this terrific help.

[ndmetherall]

Hello @swolny

I have re-run the model with a accumulation threshold of 500 instead of 1000 with the following results

Should this be ok?

I will continue work on the erodibility raster and see how it goes.

Best wishes.

[swolny]

From screenshots, the streams look fine, but what’s important is that you compare them with a real-world stream map if you can, to choose a value that creates streams that are close to reality. Since this is a simulation they will not be exact, but if needed do a few rounds of different TFA values, and see if you can get the level of tributaries pretty good. If you don’t want to run the whole SDR model to do this, you can use RouteDEM, which only does the hydrologic routing and might take less time to do a few iterations with.

No idea what happened with the erodibility raster. The layer you sent does have data in the ocean and along the coast, which makes it look like somehow you processed the datasets then clipped it in reverse, removing land and keeping water (although I was not aware that SoilGrids had soil data in the ocean at all.)

~ Stacie

[ndmetherall]

Great advice to use RouteDEM - this will save me lots of time.

From looking at the maps of streams / stream ordering of Honiara / Mataniko Catchment it looks like the TFA value of 500 may be too low since it generates more streams compared to the established river maps:

Map of wider Guadalcanal

Honiara zoomed in on Mataniko Catchment

DEM generated too many streams

So I will try a few more times in RouteDEM with higher TFA values.

In terms of the erodibility dataset I want to check I am on the right track - can you please confirm the following steps are ok? https://docs.google.com/presentation/d/1az6KH6tZzQ1vMyS7ItEAGKFnnRDgDGz9TW18aAuS1kg/edit#slide=id.ge2212007d4_0_54

1. Am I downloading the correct unit measurement datasets and rasters from Soil Grids - I am downloading SOC (dg/kg), Clay content (g/kg), Sand (g/kg) and Silt (g/kg).
   

   

2. • 0 - 5 surface depth
3. • Downloading each of the 250 m res squares one at a time (there is no way to download a larger area at one time?) - then using mosaic geoprocessing toolbox to join them all together?
     

4. I am then proceeding to do each of the following formulae to calculate K
   

5. K = fcsand * fc-si * forgc * fhisand

Eg. K = “3. Erodibility\K_equation\f_csand” * “3. Erodibility\K_equation\f_cl_si” * “3. Erodibility\K_equation\f_orgC” * “3. Erodibility\K_equation\fhisand”

Fcsand = 0.2 + 0.3 * Exp((-0.256 * “3. Erodibility\K_equation\Mosaics\SI_SD_1.tif” * (1 - “3. Erodibility\K_equation\Mosaics\SI_SL_1.tif” / 100)))

Fclay - Fcsilt = ((“3. Erodibility\K_equation\Mosaics\SI_SL_1.tif”/417)*((“3. Erodibility\K_equation\Mosaics\SI_CC_1.tif”/522) + (“3. Erodibility\K_equation\Mosaics\SI_SL_1.tif”/417))^0.3)

Forgc = 1 - (0.0256 * “3. Erodibility\K_equation\Mosaics\SI_SOC_1.tif” / (“3. Erodibility\K_equation\Mosaics\SI_SOC_1.tif” + Exp(3.72 - 2.95 * “3. Erodibility\K_equation\Mosaics\SI_SOC_1.tif”)))

Fhisand = 1 - (0.7 * (1 - “3. Erodibility\K_equation\Mosaics\SI_SD_1.tif” / 100))/((1 - “3. Erodibility\K_equation\Mosaics\SI_SD_1.tif” / 100 + Exp(-5.51 + 22.9 * (1 - “3. Erodibility\K_equation\Mosaics\SI_SD_1.tif” / 100))**

The result I tend to get looks like this which seems wrong due to the ocean values for Fcl-si: where land = 0 and sea = 1. This may be the source of the issue…?

Further details in Googledoc here: https://docs.google.com/presentation/d/1az6KH6tZzQ1vMyS7ItEAGKFnnRDgDGz9TW18aAuS1kg/edit#slide=id.ge2212007d4_0_17

[swolny]

I really can’t troubleshoot the details of your calculations, but if you’re ending up with a layer that’s all 0/1, and over the wrong areas, that seems a likely candidate. However, you are ending up with NoData over land, not 0s, so it may not be the only issue. If you’re doing these calculations in Raster Calculator, one thing to check is whether it’s creating integer values when you actually need floating point. I’ve had that happen, and it seems to help to specifically cast the problematic item(s) to float() during the calculation.

And they really don’t make it obvious if/how to download SoilGrids data for the whole globe at once, I’m sure it’s a large file but so what. One thing you can do is use WCS in a GIS session and connect to the different layers, then you can zoom and save your specific study area. Although I have really only gotten this working for SoilGrids in QGIS, not Arc for some reason (WMS works in Arc, but not WCS, and you can’t save data with WMS).

~ Stacie

[ndmetherall]

Absolutely @swolny - I understand that and would not expect anyone to have to go through all the georaster calculations.

I have been spending some time in RouteDEM and DelineateIt and have been comparing results with the existing maps:

DelineateIt (comparable but not exactly the same catchment boundaries)

RouteDEM stream features (not generating the same streams - in particulary Lungga river is not showing up).

FYI - I have tested about 6 different Flow Accumulation Values. The closest seems to be at around 10,000.

Any advice about why there might be this dissonance between the RouteDEM results and other past studies and results? Should I try to generate my streams from the ArcGIS / QGIS steps instead?

[swolny]

It is unfortunately totally typical that the streams generated from a DEM do not match reality very well. If we’re lucky, at least the major ones will match well enough, but it seems like in this case at least one (Lungga) doesn’t. It looks like the inland portion corresponds pretty well, but the outlet is in a different place along the coast. It’s probable that changing the TFA will not help with this, the TFA generally only creates more or fewer tributaries.

If the streams are not to your satisfaction, you’ve got a few options.

One is to just try a different DEM. That often helps a lot. On a recent project, I tried 4 (3 global, one national) before the 4th turned out to be satisfactory, and all created significantly different stream networks. So I recommend trying this first.

The other is to try “burning” a correct stream network into your DEM, if you have one in digital format. This usually turns out to be non-trivial, and can be very messy, and honestly I haven’t had much luck with it over the years, so am not able to give useful advice on how to do it. ArcHydro has this capability, as do other tools, or you can do it manually in GIS, trying out tutorials from the web. What I do know is that the simplistic recommendation of “subtract X large elevation value from your DEM everywhere there’s a stream” generally makes a mess for this purpose. Maybe others have had more success and can provide advice, I’d love to hear it.

There is one other work around that may (or may not) help. If you do have a digitized real world stream network, you can use it as the Drainages input to SDR, then set the threshold flow accumulation value to a very high number. That way, streams are not created from the DEM, but the model will use the Drainages input as streams, and that might work out. I’ve had mixed success with this.

Generally I greatly prefer to be working with a DEM that’s acceptable in the first place, instead of burning or doing workarounds.

Oh, one other note: If you’re using RouteDEM, make sure you’re using the MFD algorithm, since it’s the one used by SDR. If you choose D8 you won’t be getting the same streams that SDR creates.

~ Stacie

[ndmetherall]

Thanks again.

Just checking if you recommend downloading the DEM from Earth Explorer? https://earthexplorer.usgs.gov/

And which dataset exactly? I was looking at the SRTM - GTOPO30 or GTOPO30 HYDRO 1K?

However, this SRTM version should be the same as the SRTM I originally downloaded from Earth Engine.

I have also tried to download the ASTER version from the Solomon Islands Environment Data Portal.

However, when I run this through ROUTEDEM it tends to just run forever and can only generate sink filled tiff without generating the streamflow direction or stream channels and other outputs:

It has been stuck at 90.36% progress for hours:

Any advice about how I can fix this issue?

[swolny]

@ndmetherall, I don’t have a specific recommendation for DEMs, since each one seems to be better or worse, depending on the specific area in the world you’re working. For a long time, I was surprised at how badly the HydroSHEDS DEM performed in various places, until a project I’m working on now, when it turned out to be the best of 4, so it’s kind of a crapshoot I’m afraid. I’ve used SRTM and ASTER, and downloaded from USGS for various projects over the years, and none of them have been consistently the best. That said, I probably end up using SRTM most often, but not by a large margin.

I always like to try whatever DEM is provided by the government, or related provider, like the Solomon Islands Environment Data Portal, since one would hope that they’ve vetted it as being fairly correct, and it’s often trusted more by the local team (but neither of those may necessarily be the case). No idea why the pit filling is taking so long with it though. That could indicate that the DEM is particularly problematic in some way, maybe it has NoDatas scattered around or just some ornery pits. @jdouglass any idea what could cause RouteDEM to get stuck pit filling?

~ Stacie