Global River Data

For my second foray into Google Earth Engine, I decided to look at rivers, and in particular what I could do by combining the land cover data with rivers. Combining datasets opens up many new possibilities and the options quickly become overwhelming, so to keep things simple for now I aim to answer the following question: how have croplands changed near rivers over time? The question is perhaps too simple, because while it might give some indication of how productive land near rivers is, it cannot tell us whether the quality of the crops is better or worse, what kind of crops are being grown, or if the changes are compensated for elsewhere (e.g. near sources of groundwater). But it can give us an insight into changes in economic activity. If, for example, cropland near rivers disappears, what does that mean for the local farming population? what does it imply about the types of food being cultivated, or how that might affect markets? Since agriculture is often reliant on rivers, particularly in arid regions of the world, such as the River Nile in Egypt, radical changes in land use could be an indicator of a changing economic, and thus likely social situation. Of particular interest, if cropland is expanding upstream, what might that tell us about agriculture downstream, and in particular in neighbouring states?

The questions here are rather vague. We cannot hope to make economic forecasts, or geopolitical predictions, but since everything about human existence revolves around water, food and land, asking them might give insight into where we might expect to see changes as a consequence of how local economies respond to changes in both climate and neighbouring economies.

As before, this article is as much about using Google Earth Engine as the results, so I'll be discussing the technical aspects of the code as well as the data. I'll be producing maps and statistics, each product requiring its own approach.

For the cropland I will continue using the MODIS data as before, which gives us 23 years of yearly estimates of land cover. Likewise, for the geopolitical boundaries I will continue to use the geoBoundaries data. For the rivers, I chose the WWF HydroSHEDS Free Flowing Rivers Network v1 dataset (Grill et al. 2019) mainly because of the permissible licence. The only downside when compared to some other datasets available in GEE is that it was derived from a single year's data, so it will not represent changes in flow rate over time. Indeed, rivers can also change position over time, but not enough to worry me for the rather coarse statistics I will be calculating. River flow rate is seasonal, but long term trends would be worth exploring as they could represent climatic change, or change in usage upstream. I will defer this to a future test.

The key feature of the HydroSHEDS dataset for my purposes is the RIV_ORD attribute. The rivers are mapped according to their "long term average discharge" in cubic metres per second, subject to a logarithmic progression where the lowest numbers refer to the highest flow. Thus a discharge of 100,000 m³ s^-1 or more is attributed the number 1, and a discharge of between 10,000 and 100,000 m³ s^-1 is attributed the number 2 etc, down to the number 10 for discharges of less than 0.001 m³ s^-1. A map of these rivers is shown in figure 1 for RIV_ORD values 1 through 4.

We already know how to use the land cover data, so the additional steps will involve filtering to keep only cropland and masking it when it is further than a given distance from a river. The river data is contained in a feature collection where each feature consists of a LineString (a connected set of vertices tracing the path of a segment of river) with a set of properties, including the RIV_ORD attribute which we will use to filter with. These data are therefore different to the geometries and images we have been using up to now.

Since not all rivers are equal, we will only consider those above a certain discharge rate. This choice of where to draw the line is a little arbitrary, but the river networks quickly become so dense that the question of distance to a river becomes irrelevant. I tried this, and apart from the unreasonable computational load, if all rivers are mapped then even the Sahara appears to be richly served by rivers. Including these usually dry river beds is clearly absurd. As a first approximation I decided to set the cutoff at a RIV_ORD value of 4, corresponding to rivers with an average discharge of between 100 and 1,000 m³ s^-1. My litmus test was that the Murray River should be included as it counts as a significant river in south-eastern Australia, serving as a source of water, a means of transport and a state border. Many other smaller rivers in the Murray-Darling basin, such as the Goulburn, count as its tributaries and are of fundamental importance in local agriculture, though are not big enough to be included at this cutoff value. A future refinement to this work could be to descend to the next level of discharge rates. The filtering is easy enough:

// Get rivers above a certain threshold
var bigrivers = rivers.filter(ee.Filter.lessThanOrEquals("RIV_ORD", 4))

To create the mask, I made use of the convenient distance function which applies to feature collections. It generates a double precision image where each pixel's value is a measure of its distance from the nearest point, line or polygon in the feature collection. We are not yet interested in every pixel's distance, but by supplying it a cutoff distance of 100,000 m, this will define our mask, since all values beyond the cutoff are masked off.

var river_distance = bigrivers.distance(100000)  

Now we need to apply this image as a mask to the MODIS land cover data to isolate cropland within the cutoff. Because the MODIS data forms an image collection (consisting of one image per year) we need to create a function to map over each image. Here is the function definition and its evocation to the MODIS data, contained in the variable lcmode:

var mask_dist = function(image){
  var nearRiver = image.mask(river_distance)
  nearRiver = nearRiver.mask(nearRiver.where(nearRiver.neq(12), 0))
  return river_distance.mask(nearRiver)
}
var cropDistances = lcmode.map(mask_dist)

The first line in the function mask_dist takes the image and masks it using the river_distance image, such that only land cover values within the cutoff distance are preserved. The next line further masks off all land cover values that are not cropland using the where method: anything not equal to 12, which is the value identifying cropland in the MODIS LC_Type1 data, gets set to zero and is therefore masked. Now we could simply return this value and recover only cropland within 100 km of "big" rivers, but I have added an extra step, applying the mask back onto the river_distance variable. What this gives us is an image showing all cropland within 100 km of a "big" river plus its distance from the nearest river, an additional piece of information which will come in handy later.

The resultant image looks like figure 2 below.

Though this reflects the river network displayed in figure 1, let's take a moment to consider it more closely. In particular, notice that India (and indeed Europe) is quite densely covered by cropland, meanwhile the Nile's banks comprise only a thin sliver of cropland. Zooming in we find that for most of its length, there is cropland only up to a maximum distance of about 10 km either side of the river. Now comparing with India is not an apples-to-apples comparison, since India's river network is much denser and, importantly, India is also much wetter, benefitting from the South Asian monsoon. But that is the point: despite having arguably the most legendary river known to western culture, which helped propel a powerful civilisation for millennia, Egypt's food and water supply is precariously dependent on essentially this one body of fresh water, whose source is in a distant, foreign land.

The images produced in the last section can be converted into NetCDF files using, as I have, Xarray in Python, and exploited however one likes. While they might be interesting to look at, it is more useful to generate statistics for some quantitative analysis. This is best done within Google Earth Engine. To begin, we do much the same as for the script for the yearly maps by filtering the cropland within a threshold distance of the rivers with the distance method. Then, just as for the land cover statistics, we use reduceRegions to calculate statistics for a given region. This time I didn't just want a simple sum, I wanted to take advantage of the distance function attributed to the cropland pixels. It would therefore be nice to perform a sum for distinct distance intervals from the rivers. The ee.Reducer class includes histogram functions, of which fixedHistogram is the most appropriate for our needs. By predefining the intervals relative to the cutoff distance thusly,

var minDistance = 0
var maxDistance = 100000
var sliceWidth = 10000
var slices = (maxDistance-minDistance)/sliceWidth

we can construct a histogram through reduceRegions thusly,

var genHist = function(image){
  var histograms = image.reduceRegions({
      reducer: ee.Reducer.fixedHistogram({
      min: minDistance,
      max: maxDistance,
      steps: slices,
      cumulative: true,
    }).splitWeights(),
    collection: geolabFiltered,
    scale: 1000
  })
  return histograms
}

The splitWeights() method is used to weight the sums by the pixel area, and it requires some explanation, partly because it was unnecessary in the land cover sums from before (they were simply sums of pixels multiplied by area), but especially since the GEE Reference manual's description of it is so unhelpful:

Returns a Reducer with the same outputs as the given Reducer, but with each weighted input replaced by two unweighted inputs. (Google for Developers n.d.--a)

The Guide is slightly better, stating that

A reducer modified by splitWeights() takes two inputs, where the second input is the weight (Google for Developers n.d.--b)

But how those inputs should be arranged might not be clear to someone not used to GEE's paradigm, though it is simple enough: the weights, in this case the area per pixel, should be the second band in the image. There is no way to specify the band explicitly as an argument to splitWeights so this requires us to manually place the weight in that second position, if it isn't (as in my case) already. This too is easy enough, one simply selects the bands in the desired order. The generation of the pixel areas as an extra band is as before

var add_country_band = function(feature){
  return ee.Image.pixelArea().addBands(feature)
} 
var riverDistArea = cropRiverDist.map(add_country_band)

And we select the bands with "area" second with

riverDistArea = riverDistArea.select(["distance", "area"])

We can now map the genHist function defined above onto riverDistArea and splitWeights() will do its job,

var hist = riverDistArea.map(genHist)

The end result is a histogram of land cover according to distance intervals from rivers, per country. I chose to do cumulative sums, and set the scale to 1 km. Here is an extract of the table with the cumulative areas as columns for each country, for every year of the MODIS data, where I've only shown the columns for distances up to 50 km from a river. Units are in square metres.

Of course these statistics are most interesting when viewed as a function of time. Seeing as we've discussed it, let's look at Egypt's statistics for the full period.

I've converted the units from metres (for the distance labels) and square metres to kilometres and square kilometres. Each column shows the surface area of cropland for the 10 km interval starting at the column's heading, so the first column shows the area up to 10 km from rivers, the second between 10 and 20 km and so on. The "total" column is the sum of all other columns, and the "countrywide" column shows the total land cover for all of Egypt, whose calculation was discussed in the previous article. We can more easily read this as a plot where we subtract the 2001 value from each column.

This gives us an indication of each distance interval's total gain and rate of change. It shows that over the 22 years of the MODIS data, cropland increased in surface area in every one of those intervals, and that indeed most of the cropland is within 100 km of the Nile or its Delta. Some caution should be exercised in interpreting exactly how much lies outside the 100 km limit, since recall that we have only measured from large rivers, and it is entirely possible that many small streams are used for crop irrigation. It is also possible (and indeed a quick visual verification bears this out) that the cropland outside the 100 km limit is still close to it. By consulting a map, we find that most of the cropland is in the Delta. Either way, the proportion outside the limit (comparable to the last few intervals) is small enough that it won't have a major effect on the conclusions we derive. I'll come back to the Delta, because although much of it is criss-crossed with canals, these canals do not all show up in the river data.

The greatest increase in cropland is in the 10–20 km interval. This reflects two things: that land close to the river is more valuable cropland (easier to irrigate), but also that the land closest to the river is probably saturated. Indeed, from the table we see that the number of square kilometres of cropland in the interval closest to the river is well over double that of this second interval even at the end of the period, and that the second interval comprises almost a third of the total increase up to 100 km. It is interesting to find that some intervals actually lose substantial amounts of cropland at various periods. This could be because the land is temporarily devoid of crops or for some reason the signature is unclear and fools the data processing algorithm, or alternatively the land is put to use for something else and later on new cropland is opened up. Without looking at the satellite images in detail this is not possible to ascertain.

Let us leave aside for a moment the question of distance to major water ways, whose application is, admittedly, limited. There are many reasons to release more land for agriculture: economic incentives to produce more agricultural products, such as for export, or to compensate for possibly less reliable growing seasons in a changing climate. There are also several factors that may contribute to the feasibility of generating more cropland, such as improved water management, increased fresh water transport by the Nile, more access to labour or more efficient farming techniques. But next to these incentives we must give centre stage to Egypt's positively exploding population—from the beginning to the end of the MODIS data period it increased from around 70 million to approaching 120 million, and it is currently climbing by about 2 million people per year (World Population Review n.d.). The relatively modest increase in cropland, about 8% from 2001 to 2023, on its own would not be sufficient to feed these extra people, so either the farming would have to have become much more efficient, or the import–export relationship would have to have changed substantially; or perhaps a combination of the two. It seems unlikely to me that farming efficiency has increased that much in just 22 years, so the economic and social implications of this increase are dramatic. With the population set to increase still further, Egypt is certainly facing some challenges. But if the current landscape is anything to go by, with all the cropland (and population) thus far being concentrated so close to the Nile, and moreover dependent on it as the only major supply of fresh water, it would appear that expansion of food production will be among the greatest of these challenges.