Hack 76. Explore the Effects of Global Warming
See what sea-level change means for the world's coastlines and the people who live near them .
If present trends in greenhouse gas production, deforestation, and average surface temperature continue, it seems likely that Earth's ice caps will continue to melt, dumping untold
6.14.1. Importing the Elevation Data
The first part of this questionwhere will the new coastlines bedepends on where the current coastlines are. We can start by getting a one-kilometer resolution digital elevation model (DEM) from the GLOBE Project and importing it into GRASS. The GLOBE project's site is http://www.ngdc.noaa.gov/seg/topo/globe.shtml, but you can get the data directly from http://www.ngdc.noaa.gov/seg/topo/globeget.shtml. Click the link to "Select Your Own Area," and you'll be taken to a form where you can specify the details of the DEM you want. Select the Custom... region from the drop-down box on the left and then enter the bounding box for the DEM. We selected -12 degrees for the western edge, 30 degrees for the
The DEM will be in ESRI Binary Interleaved, or BIL format , which can fortuitously be read by GRASS's r.in.gdal tool. Make a new directory, and unpack the tarball inside it:
GRASS 5.7.0:~ > mkdir climate GRASS 5.7.0:~ > cd climate GRASS 5.7.0:~/climate > tar xvfz ~/Desktop/europe_dem.tgz europe_dem.bil europe_dem.hdr europe_dem.clr GRASS 5.7.0:~/climate > r.in.gdal -o in=europe_dem.bil out=dem 100% CREATING SUPPORT FILES FOR dem
The tarball contains three files: a
file that contains the actual elevation data, a
file that contains the metadata, and a
file that contains a default
We can assign a simple green-yellow-red color map to the dem layer, set the current region to that of the raster layer, and then display it on an X monitor:
GRASS 5.7.0:~/climate > r.colors dem color=gyr Color table for [dem] set to gyr GRASS 5.7.0:~/climate > g.region rast=dem GRASS 5.7.0:~/climate > d.mon start=x0 using default visual which is TrueColor ncolors: 16777216 Graphics driver [x0] started GRASS 5.7.0:~/climate > d.rast dem bg=blue 100%
Since the DEM layer import sets existing ocean
Figure 6-47. A
Now, if we want to visualize which
-103 blue 10 blue 11 green 1000 yellow 4570 red
The new set of color map rules says "
GRASS 5.7.0:~/climate > r.colors dem color=rules <dem.rules Color table for [dem] set to rules GRASS 5.7.0:~/climate > d.redraw 100%
The first thing we notice is that the Netherlands almost disappear completely! Additionally, major cities like London, Copenhagen, and Venice fare rather poorly. This visualization practically begs our other question: if the sea level were to rise by 10 meters, how many people might be displaced? To address this question, we'll have to redo our analysis to
Our second approach will use
to identify which parts of our elevation model lie below 10 meters, and hence are at risk of being submerged. Raster algebra is simply any process that derives each cell of a new output layer by applying a
GRASS 5.7.0:~/climate > r.mapcalc 'submerged = if(dem <= 10, 1, null())' 100% GRASS 5.7.0:~/climate > r.colors submerged color=rules Enter rules, "end" when done, "help" if you need it. Data range is 1 to 1 > 1 blue > end Color table for [submerged] set to rules
Our call to r.mapcalc creates a new layer called submerged , which contains a value of 1 everywhere dem contains a value less than or equal to 10, and a null value everywhere else. This gives us a raster layer that consists solely of the land areas that would be submerged, which we then color blue. If we reset the colors of our dem layer, we can display the original coastlines subsequently overlaid with the new coastline, giving us a map that looks schematically similar to the map we made in Figure 6-49:
GRASS 5.7.0:~/climate > r.colors dem color=gyr Color table for [dem] set to gyr GRASS 5.7.0:~/climate > d.erase GRASS 5.7.0:~/climate > d.rast dem bg=blue 100% GRASS 5.7.0:~/climate > d.rast -o submerged 100%
The -o option to d.rast causes it to filter out the null values in submerged , rather than erasing the dem layer underneath, resulting in an overlay effect. If we want to highlight, rather than blend in, the submerged areas, we can assign a different color to that layer and redraw. Figure 6-50 shows the highlighted areas overlaid on the original map.
GRASS 5.7.0:~/climate > r.colors submerged color=rules Enter rules, "end" when done, "help" if you need it. Data range is 1 to 1 > 1 red > end Color table for [submerged] set to rules GRASS 5.7.0:~/climate > d.redraw 100%
All of the potential elaborations of this idea, such as using r.mapcalc to combine the original and submerged land areas into a single layer or using v.in.ascii (or s.in.ascii in GRASS 5.3 and earlier) to import and overlay a map of major European cities, are left as exercises for the reader, so that we can get back to our second question.
At last, we return to our second question: how many people in Western Europe might be affected by a 10-meter rise in global sea level? To answer this, we will need to import a layer containing population data, which we can combine with our elevation layers to begin making estimates. The "Gridded Population Map of the World,"
GRASS 5.7.0:~/climate > unzip v2eup95agi.zip Archive: v2eup95agi.zip inflating: eup95agi.bil inflating: eup95agi.blw inflating: eup95agi.hdr extracting: eup95agi.stx inflating: readme.txt
However, in practice, the GDAL driver has no way to distinguish between 32-bit floating point values and the 32-bit integers actually stored in the BIL file, due to a bug in the BIL format specification. There are two ways around this. One is to use r.in.bin , but that means calculating the extents of the file by hand. The other way is to use GDAL's "virtual raster" driver to explicitly specify the file structure of the BIL file by means of a separate XML file. We'll take the latter approach. If you run gdalinfo on the data set, you can get the extents and spatial resolution (i.e., degrees per pixel) of the raster layer. We'll use the -mm option to have GDAL count the minimum and maximum values, just to verify that it has trouble with this data set:
GRASS 5.7.0:~/gis/europe > gdalinfo -mm eup95agi.bil Driver: EHdr/ESRI .hdr Labelled Size is 1680, 1320 Coordinate System is `' Origin = (-25.000000,85.000000) Pixel Size = (0.04166667,-0.04166667) Corner Coordinates: Upper Left (-25.0000000, 85.0000000) Lower Left (-25.0000000, 30.0000000) Upper Right (45.0000000, 85.0000000) Lower Right (45.0000000, 30.0000000) Center (10.0000000, 57.5000000) Band 1 Block=1680x1 Type=Float32, ColorInterp=Undefined Computed Min/Max=0.000,0.000
As you can see, GDAL gets the minimum and maximum values wrong, as well as the data type ( Float32 instead of Int32 ). Using the geographic origin, and the layer height, width, and pixel size indicated by gdalinfo , we can create a GDAL virtual-raster description that explicitly sets the cell data type to 32-bit integer. Fire up a text editor and dump the following to a file called eup95agi.vrt :
<VRTDataset rasterXSize=" 1680 " rasterYSize=" 1320 "> <GeoTransform> -25 , 0.0416666667 ,0, 85 ,0, -0.0416666667 </GeoTransform> <VRTRasterBand dataType=" Int32 " band="1" subClass="VRTRawRasterBand"> <SourceFilename relativetoVRT="1"> eup95agi.bil </SourceFilename> <ByteOrder> MSB </ByteOrder> </VRTRasterBand> </VRTDataset>
As you can see, we've
Now you can treat this XML file as if it were an ordinary GDAL layer, by running the following in GRASS:
GRASS 5.7.0:~/climate > r.in.gdal -o in=eup95agi.vrt out=population 100% CREATING SUPPORT FILES FOR population
In our new population layer, each cell has a value equal to the estimated number of people that lived in the area represented by that cell in 1995. Unlike the elevation data, the population data doesn't differentiate between water areas and places that are simply uninhabited, which would appear as an undifferentiated mess if you ran d.rast population at this point. (Try it and see, if you like. d.what.rast can be used to query individual raster cells with your mouse.)
To address this, we'll
GRASS 5.7.0:~/gis/europe > r.null population setnull=0 Writing new data for [population]... 100% GRASS 5.7.0:~/climate > r.colors population color=rules Enter rules, "end" when done, "help" if you need it. Data range is 0 to 626281 > 0 green > 1000 yellow > 10000 red > 700000 purple > end Color table for [population] set to rules GRASS 5.7.0:~/climate > d.erase GRASS 5.7.0:~/climate > d.rast population bg=blue 100%
Once the custom color map is applied, we get a very nice looking population map of Europe, as shown in Figure 6-51.
Now, back to our second question, which was, "How many people in Western Europe will be displaced from their
As a check, let's run r.sum on our population layer and see if it produces a reasonable value. In order to make this work, we'll need to make sure that the resolution of the current GRASS region matches that of the population layer, or some values might end up being counted multiple times, which will throw the sum off. Running r.info on the population layer yields the following metadata:
GRASS 5.7.0:~/climate > r.info population ... N: 85N S: 30N Res: 0:02:30 E: 45E W: 25W Res: 0:02:30 ...
This tells us that the spatial resolution of the image is 2.5 minutes per pixel. We can feed this value directly into g.region and then run r.sum on the population layer:
GRASS 5.7.0:~/gis/europe > g.region res=0:02:30 GRASS 5.7.0:~/gis/europe > r.sum population Reading pop... 100% SUM = 604283674.000000
The value we get back, around 604 million, seems about right for the area depicted on our graphics monitor. Now let's create a displaced layer that contains the population values for all the areas in the submerged layer. In order to do this, the submerged layer and the population layer have to have the same resolution, or we'll be comparing apples to oranges. We can use r.resample to down-sample our submerged layer to the resolution of the population layer, and then use r.mapcalc to synthesize the displaced layer:
GRASS 5.7.0:~/gis/europe > r.resample in=submerged out=submerged2 percent complete: 100% Creating support files for submerged2 creating new cats file... GRASS 5.7.0:~/gis/europe > r.mapcalc 'displaced = if(isnull(submerged), null(), population)' 100%
Finally, we can use r.sum to estimate how many people live in the potentially submerged areas:
GRASS 5.7.0:~/gis/europe > r.sum displaced Reading displaced... 100% SUM = 41109241.000000
Our answer is that about 41 million people will be displaced from their homes in Western and Central Europe, should the current sea levels rise by 10 meters. That's more than one in every 15 people! If we perform the same analysis for the East