Watersheds as art

I was delighted by the look of this watershed. Can you figure out where it is?

Just a few kilometers west of Bolivia’s capital, Sucre, lies the Maragua Crater, the site of a long-ago asteroid impact.

All of the top search results are about hiking out to visit it. There are 2,000 year old cave paintings, fossil dinosaur footprints, and incredible multicolor rock formations. Wow, new place to add to the bucket list!

Photo by Cody Hinchliff, flickr.

Web hosting expenses have gone up

I recently became aware that the Global Watersheds app occasionally goes down when it’s under heavy load. At its peak, the app is delivering a watershed every 1.2 seconds. You all are loving the site to death! 🙂

I’ve upgraded my hosting plan to 2GB of RAM. Let’s see how much that helps. It costs a few hundred dollars per year, which is a lot for a hobby project by a full-time student. If you’ve enjoyed using the app, or it’s been valuable in your work, would you please support it by sending a few dollars? Thank you! 🙏 💕

Sharing watersheds with the world

The Global Watersheds web app has a nice feature that it seems only a few people have discovered. You can create a permalink to an interactive map of your watershed (or flowpath). Then you can embed the map on a web page, send it to a friend, or post it on your favorite social media site. Just look for these buttons at the bottom left after you have created a new watershed:

Since the app was launched in October 2022, users have created 770 shared watersheds. Here is a map and a table showing all of them. In that same time, users have made over 237,000 watersheds and over 36,000 flow paths.

Happy exploring!

Decreasing the white space between Matlab subplots — the easy way

Matlab is great for producing high-quality graphs and plots for all kinds of science and engineering problems. Plots can exported to many formats, and inserted directly into reports or journal articles. 

Matlab is also great when you want to create multiple plots together on a single figure. The standard way to do this is with the subplot() command. Indeed, you can create attractive figures with a grid of plots. The following code creates a 3 x 3 set of 9 plots:

figure('units','normalized','outerposition',[0 0 1 1]);
for i = 1:9
   subplot(3, 3, i)
   plot(rand(1, 20));
end

A common complaint about the subplot function is that it leaves too much whitespace in between the plots. This is especially noticeable when you are making plots without tick labels on the x-axis or y-axis: 

tiledlayout()

In 2019, Matlab introduced the tiledlayout() function as an alternative to subplot(), which gives you more control over the amount of space between the individual plots. Here’s an example of it’s use:

figure('units', 'normalized', 'outerposition', [0 0 1 1]);
t = tiledlayout(3, 3, "TileSpacing", "tight");
for i = 1:9
   nexttile
   plot(rand(1, 20));
end

The ‘TileSpacing’ parameter can take the values ‘tight’, ‘none’, ‘compact’, or ‘loose’. This gif shows the difference: 

tight_subplot()

The tiledlayout function works well, most of the time. If you need even more control over the layout of your subplots, head over to the Matlab File Exchange, and download the function tight_subplot(). There are many similar contributions, but this one appears to be the most popular, with 47K downloads as of early 2023. It works perfectly and is simple to use once you get the hang of it.

With the tight_subplot() function, you can specify exactly how much space you want between plots. 

You can separately adjust the vertical spacing and horizontal spacing between subplots. 

To do this, adjust the parameters gap_h and gap_w for the gap distance between plots in the horizontal and vertical directions. 

[ha, pos] = tight_subplot(rows, cols, [gap_h, gap_w]);

The function returns two variables. The first is an array of handles of the axes objects. You have to activate the axis in order to make a plot. For example:

[ha, pos] = tight_subplot(3, 3, [0.05, 0.05]);
for i = 1:9
    axes(ha(i));     % Don't forget to 'activate' the axis. 
    plot(rand(10, 1); 
end

The values for gap_h and gap_v are “normalized units” from 0 to 1, where 1 is the full width or height of the computer screen is 1. So setting gap_h to 0.01 means the horizontal gap between the plots will be 1% of the screen width. Here is an illustration of changing gap_h:

Here is the effect of modifying gap_w: 

And here’s an example where I used the tight_subplot function to create a set of maps of global evapotranspiration over 9 consecutive months, using data from a meteorological reanalysis model called ERA5. This works well because we don’t need axis labels. I also used the trick of creating a single colorbar and customizing its position on the figure.

How to calculate the average precipitation over a watershed with gridded data

In this post, I’ll give a short tutorial for how to calculate the average precipitation over a watershed or river basin. This is a common task in hydrology and environmental science.  

Historically, hydrologists used data from rain gages, which report precipitation at a single location. Yet, everyone knows that rainfall varies a lot from one place to another. So for a large watershed, you might have gathered data from several gages. Hydrologists have developed several clever methods for averaging point observations from multiple gages. 

Today, gridded climate data are widely available. For precipitation, gridded data can give you much more information about how rainfall varies with geography. Also, they can give you information for remote or sparsely populated regions where rain gauges are scarce.

If you’re in the US, you might be using PRISM, which is based on a sophisticated interpolation of data from hundreds of gages. Or you could be using a global dataset based on satellite remote sensing, for example CMORPH. If you are looking back at history and need long records, you might choose output from a reanalysis model like NCEP or ERA5

Below, I show how to calculate the spatial mean of these data. I use precipitation as an example, but the same methods will work with any kind of gridded environmental data – evaporation, temperature, land use, vegetative cover, etc. I’m also talking about watersheds, but the same methods could be used to get the average over a city, a province, the boundaries of a bioregion, etc. 

If you only need to do this calculation once, you can use GIS to calculate a zonal average. I’ll show how to do it with the free software QGIS. If you need to do this calculation many times (i.e. with daily precipitation), you will want to write code to automate this. I’ll show how to do that in a future post. 

Example Application: Flooding on the Winooski River

Here, we’ll estimate the amount of rain that fell over the Winooski River watershed in Vermont on July 11, 2023, a day where there was major flooding. Here are the steps:

1. Get PRISM precipitation

Go to https://prism.oregonstate.edu/. There are lots of different options. I downloaded provisional daily precipitation for July 11, 2023. Here’s what it looked like: 

Unzip the files to a convenient location. The data are in BIL format. This is an old ESRI format for aerial photos and remote sensing data, but it shouldn’t pose a problem if you have a full installation of QGIS. 

2. Get your watershed boundary

Go to the global watersheds web app at https://mghydro.com/watesheds

I panned and zoomed until I found where the Winooski drains to Lake Champlain. 

Under Options, check the box for “Make results downloadable.” 

Click on the map then click “Delineate!” button in the map popup. Or you can click “Enter coordinates” and enter 44.53, -73.27. It should look something like this:

If the results don’t look right, click in a slightly different location and try again. 

On the left of the page, scroll down. Under Downloads, click Watershed Boundary. Click the button to download the watershed boundary. I recommed choosing a GeoPackage, but the other formats will work fine too. 

3. Create a map in QGIS

Open QGIS, and create a new project. 

Add the watershed: Select Layer > Add layer > Add vector layer, then choose the watershed layer.

Add the precipitation layer: Select Layer > Add layer > Add raster layer, then choose the PRISM precip. layer. Choose the .bil file.

Here, I adjusted the Symbology of the layers to make them look nice.

Use the “Identify features” tool to check a few values of the precip. We can see a pixel in the center of the watershed where the precip was 95 mm on July 11. That is about 3.7 inches. A lot of rain in 24 hours! 

4. Calculate the basin average

Open the Toolbox It looks like a little gear in the menu bar, or choose Processing > Toolbox. 

Search for the tool Zonal Statistics, and double click it to open. We have to make the right selections in the window that pops up:

Under Input layer, select the watershed vector layer. 

Under Raster layer, choose the PRISM precipitation raster. 

Under Raster band, keep the default, Band 1. (This raster only has one band. Sometimes a raster will have multiple bands. For example, an image will have separate bands for Red, Green, and Blue.)

Under Statistics to calculate, make sure Mean is included. 

Near the bottom, you can keep [Create temporary layer], or you can choose to save the results. Your choices are a variety of geodata layers. Since the results will be a table (not geodata), I  chose .csv, a comma-delimited text file. 

Click Run. 

In a moment, you should see a new table appear in your map’s Table of Contents. Right click on it and choose Open Table. 

Note that the table has only one row. That is because our input vector file only had a single feature. 

In the table,-mean is 56. That means the watershed received an average of 56 mm of precipitation that day. The field _count has a value of 176. That means that QGIS averaged the value of 176 pixels that itersect our watershed. 

Next Steps

That’s it! Now you know how to calculate the average precipitation over a watershed. This kind of calculation is extremely important in many areas of science and engineering. It’s useful for analyzing floods and droughts, in water budget studies, etc. 

The approach we used required a lot of clicking. If you need to do it over and over, you can write some code to automate the calculation. Let me know if you’re interested in seeing this in a future post. 

Important note: This method, using the zonal average in GIS, works well when your watershed is small. That is because the pixels all have roughly the same area. If you are dealing with a large watershed, the results will not be accurate, because the area of the pixels varies a lot with latitude. This illustration shows how grid cells get much smaller toward the poles.

For larger watersheds, you should calculate a weighted average that accounts for the varying area of the pixels. This means you’ll have to write some code to do it.

Mapping the Cuyahoga River’s Watershed

The Cuyahoga River in Ohio is famous to anyone who has studied the history of the environmental movement in the United States. Outside Magazine published a great article about the river in 2020, recounting its pivotal role in the history of the Clean Water Act, and describing ongoing efforts to clean up the river.

Here is a map of the Cuyahoga watershed using data from the US Geological Survey. I recently added this as a new data source option in the Global Watersheds mapper.

Improved Flowpath Tracing

Last month, I posted a poll on the Global Watersheds page. The app is getting around 200 users per day, and I was curious to see who is using it and why.

I found out that most of you are just as interested in tracing downstream flowpaths as you are in delineating watersheds. To be honest, I added this feature to the app as a bit of an afterthought.

More detailed flowpaths

Since folks are so interested in downstream flowpaths, I decided to add a couple of new features. You can now use higher-precision mode. (This is enabled by default if you are zoomed in far enough or when you manually enter coordinates).

Previously, the flowpath started at the closest river or stream. Now, the flowpath starts right near the point where you clicked. (Well, starting at the nearest pixel with a downstream neighbor.) The app will trace the flow from the point of origin, going downhill over land until it encounters a stream or river. From there, it flows to the sea or to an inland sink.

Share flowpaths

You can now share or embed your downstream flowpaths. After you’ve created a new flowpath, just click one of the sharing buttons at the left bottom. The link will point to a custom page with your flowpath. To embed a map widget, click the embed button </>, and copy and paste the code into a web page or blog post.

Tracing contaminant flow

Here’s an example of both new features in action. On February 2, 2023, a freight train derailed in East Palestine, Ohio. According to the US Environmental Protection Agency, vinyl chloride and other toxic chemicals were “released to the air, surface soils, and surface waters.”

Based on the EPA reports, the spill occurred at latitude: 40.836° N, longitude: -80.522° W. We can use the Global Watersheds app to find out where water flows downstream from this location. This can help us begin to understand where these chemicals can be transported by water, and what people and ecosystems could be at risk.

As always, if you encounter any bugs using these new features, send me a message!

Some new map options in the Global Watersheds web app

On the Global Watersheds app, I noticed that the OpenTopo map layers were not loading for me at all. I think I may have gone over some kind of limit.

As a workaround, I’ve added a whole bunch of new basemap options. A couple of them only work in certain regions. For example, GeoPortail France only works at higher zoom levels in France, and the USGS topographic maps can only zoom in to the US.

Some of the new basemap options have limits on the number of map tiles they will provide for free. So if one doesn’t work, try a different one. And if you have any feedback, let me know! 🙂

Demo: Use the Global Watersheds API and Python to automatically delineate watersheds

It seems that a few people are discovering the mghydro watershed delineation API. Last Monday, the API handled requests for 9,679 watersheds over the course of a few hours. (Someone is very interested in the watersheds of France!)

One way to use the API is with GIS software, to fetch watershed boundaries and river centerlines from the web app and display them in your desktop GIS. Another use is to automate your requests, rather than using the online map interface. This could be useful if you want to delineate hundreds (or thousands!) of watersheds. However, remember that the results of automated delineation are very often wrong, so you should visually check each one.

Here is a demo in Python of how you could use the API to download watersheds for multiple outlet points. Click on the filename at the bottom to go to Github, where you can download the Jupyter notebook.

The Genesee River Watershed

I grew up a stone’s throw from the Genesee River, in Rochester, New York, but never knew that the river began in far-off Pennsylvania. Exploring watersheds can change the way you view the world.

According to family lore, my ancestors (around the 1820s or 30s?) were headed to the western frontier, but when they got as far as Rochester, the women declared they were hot, dusty, and tired, and this little town with its flour mills and waterfalls suited them just fine.

Full screen: https://mghydro.com/watersheds/shared/F15986.html

As the Genesee flows north through central New York State, it carves out the magnificent Grand Canyon of the East, which you can visit at Letchworth State Park. This was a popular day trip for my family when I was young, and must have influenced my lifelong love of rivers!

The Genesee River gorge in Letchworth State Park, with the historic 1875 Erie Railroad bridge crossing upstream of Upper Falls. Photo by James St. John.