Module 3: Extreme weather
Script 02: Future Climate Projections
Max Hofmann
Welcome to the second script of Module 3! Here, we want to investigate the amount of heat events forecasted for the Caucasus region for the future. At the example of future climate projections, you will here learn how to:
- work with NetCDF files
- identify extreme weather events
- summarize the occurrence of extreme weather events in space and time
Please read the following code instructions carefully, try to understand the code that follows each instruction, execute it and see what happens. Do not hesitate to change the code or write your own code, and execute it as well!
1. Import NetCDF files
We start off again with defining our working directory and loading all the packages that we will need:
dir <- "C:/Users/Max/Desktop/IAMO_neu/eLearning/Module_3_Extreme_weather/M03_Part02_Future-Climate-Projections/"
library(reshape2)
library(dplyr)
library(ggplot2)
library(ncdf4)
library(raster)
library(rgdal)In the folder “data”, there are 16 NetCDF files that contain daily future estimates for maximum air temperature (tasmax) in the Caucasus region with a resolution of 0.5 degrees. There is one file for each RCP (RCP4.5 and RCP8.5), future time period (2051-2060 and 2081-2090), and climate model (GFDL-ESM2M, HadGEM2-ES, IPSL-CM5A-LR and MIROC5). NetCDF files can be opened with the function nc_open() from the package ncdf4. We import the file for MIROC5, RCP8.5, 2081-2090. For simplicity, we call this object “miroc” here:
miroc <- nc_open(paste0(dir, "data/tasmax_day_MIROC5_rcp85_r1i1p1_EWEMBI_20810101-20901231_lat38.0to45.0lon39.0to52.0.nc4"))Let’s have a look at the resulting object which is of type “ncdf4”. We see that it has 1 variable, which is air temperature in degrees Kelvin, called “tasmax”. Each measurement of “tasmax” has 3 associated values pertaining to the three dimensions: longitude (lon), latitude (lat), and time (time) expressed as days since January 1, 2006. We also find in “miroc” some metadata, e.g. information about the publishing institution and the project in the frame of which the dataset was created:
miroc## File C:/Users/Max/Desktop/IAMO_neu/eLearning/Module_3_Extreme_weather/M03_Part02_Future-Climate-Projections/data/tasmax_day_MIROC5_rcp85_r1i1p1_EWEMBI_20810101-20901231_lat38.0to45.0lon39.0to52.0.nc4 (NC_FORMAT_NETCDF4_CLASSIC):
##
## 1 variables (excluding dimension variables):
## float tasmax[lon,lat,time] (Chunking: [26,14,1]) (Compression: shuffle,level 1)
## _FillValue: 1.00000002004088e+20
## missing_value: 1.00000002004088e+20
## standard_name: air_temperature
## long_name: Daily Maximum Near-Surface Air Temperature
## units: K
##
## 3 dimensions:
## lat Size:14
## axis: Y
## standard_name: latitude
## long_name: latitude
## units: degrees_north
## lon Size:26
## axis: X
## standard_name: longitude
## long_name: longitude
## units: degrees_east
## time Size:3652 *** is unlimited ***
## standard_name: time
## units: days since 2006-1-1 00:00:00
## calendar: proleptic_gregorian
## axis: T
## long_name: time
##
## 6 global attributes:
## Conventions: CF-1.6
## title: CMIP5 output of tasmax_day_MIROC5_rcp85_r1i1p1, first-order conservatively interpolated to a 0.5 degree regular grid and bias-corrected using EWEMBI data from 1979 to 2013
## institution: Potsdam Institute for Climate Impact Research, Research Domain Climate Impacts and Vulnerability, Potsdam, Germany
## project: Inter-Sectoral Impact Model Intercomparison Project phase 2b (ISIMIP2b)
## contact: ISIMIP Coordination Team, Potsdam Institute for Climate Impact Research (info@isimip.org)
## references: See ISIMIP2b bias-correction fact sheet (http://www.isimip.org/gettingstarted)class(miroc)## [1] "ncdf4"The object “miroc” is not ready to work with yet. To work with a variable or dimension contained in “miroc” that we are interested in, we have to first import it with the function ncvar_get(). Let us import the temperature variable “tasmax” and name the corresponding object “miroc_temp”. This object is an array that consists of 26 rows, 14 columns, and 3652 layers:
miroc_temp <- ncvar_get(miroc, "tasmax")
class(miroc_temp)## [1] "array"dim(miroc_temp)## [1] 26 14 3652There are 26 rows and 14 columns in “miroc_temp” because our dataset covers the extent from 38° to 45° North and 39° to 52° East, with a spatial resolution of 0.5 degrees. Accordingly, “lon” and “lat” are vectors of 26 and 14 values that contain the center coordinates of the respective cells:
miroc_lon <- ncvar_get(miroc, "lon")
miroc_lon## [1] 39.25 39.75 40.25 40.75 41.25 41.75 42.25 42.75 43.25 43.75 44.25 44.75
## [13] 45.25 45.75 46.25 46.75 47.25 47.75 48.25 48.75 49.25 49.75 50.25 50.75
## [25] 51.25 51.75miroc_lat <- ncvar_get(miroc, "lat")
miroc_lat## [1] 44.75 44.25 43.75 43.25 42.75 42.25 41.75 41.25 40.75 40.25 39.75 39.25
## [13] 38.75 38.25After importing all variables, you should close the NetCDF file with the function nc_close():
nc_close(miroc) We could now proceed working with the temperature data in “miroc_temp”, which is in the format of an array. However, NetCDF files can also be imported as RasterStacks:
miroc_stack <- stack(paste0(dir, "data/tasmax_day_MIROC5_rcp85_r1i1p1_EWEMBI_20810101-20901231_lat38.0to45.0lon39.0to52.0.nc4"))Why do we then actually need NetCDF files? TIFF files are somewhat limited, whereas NetCDFs can attribute more than 3 dimensions to your data, and contain several variables at once. They can also contain more metadata than a TIFF file, as you have seen above. However, the NetCDFs from ISIMIP are not very complex and for our current purposes, it is sufficient to import them as RasterStacks.
2. Identify heat events
Our goal here is to identify future heat events in Armenia. We define such an event as one raster cell exceeding a certain temperature threshold on one day. When two raster cells exceed the threshold during the same day, that counts as two heat events. When one raster cell exceeds the threshold during two days, that also counts as two heat events. In the end, we calculate the amount of heat events for large periods of time to investigate whether there is a temporal trend in the modeled future occurrence of heat events.
Let us have a closer look at “miroc_stack”. The dimensions and the extent match our observations from the ncdf4 object “miroc” that we inspected abvoe. However, in “miroc_stack”, the names are already actual dates, which is helpful for us:
miroc_stack## class : RasterStack
## dimensions : 14, 26, 364, 3652 (nrow, ncol, ncell, nlayers)
## resolution : 0.5, 0.5 (x, y)
## extent : 39, 52, 38, 45 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2081.01.01, X2081.01.02, X2081.01.03, X2081.01.04, X2081.01.05, X2081.01.06, X2081.01.07, X2081.01.08, X2081.01.09, X2081.01.10, X2081.01.11, X2081.01.12, X2081.01.13, X2081.01.14, X2081.01.15, ...The stack contains 3652 layers, which is equivalent to 10 years of daily data. The first day in the stack is January 1, 2081, and the last one is December 31, 2090:
names(miroc_stack)[1]## [1] "X2081.01.01"names(miroc_stack)[3652]## [1] "X2090.12.31"For simplicity, We select the first layer of “miroc_stack” and work with this one for a while:
miroc_sample <- miroc_stack[[1]]We import the Armenia shapefile and plot it on top of “miroc_sample” to get an idea what area is covered by the stack:
armenia <- readOGR(paste0(dir, "shapefiles/gadm41_ARM_0.shp"), verbose=F)
plot(miroc_sample)
plot(armenia, add=T)
You can see from the legend in the map above that the temperature measurements are in degrees Kelvin. We can easily change that with the following formula:
miroc_sample_C <- miroc_sample - 273.15
miroc_stack_C <- miroc_stack - 273.15
plot(miroc_sample_C)
plot(armenia, add=T)
In module 2, you have already learned how to classify a raster with logical expressions. Let us create a new raster that is “0” where “miroc_sample_C” is below or equal to 10 °C, and “1” where it is above 10 °C. We see that Armenia intersects with two of such cells on January 1, 2081:
miroc_sample_C_heat <- miroc_sample_C
miroc_sample_C_heat[miroc_sample_C_heat <= 10] <- 0
miroc_sample_C_heat[miroc_sample_C_heat > 10] <- 1
plot(miroc_sample_C_heat)
plot(armenia, add=T)
Let us have a look how we can automatically count these cells!
We first have to cut out Armenia from the whole region. When we simply apply the mask() function with the shapefile, the problem is that some cells that are only marginally covered by Armenia get lost:
plot(mask(miroc_sample_C_heat, armenia))
plot(armenia, add=T)
We can circumvent this by converting the shapefile “armenia” to a raster. The resulting raster “armenia_rastered” contains values for all cells intersecting with the shapefile, and the cell values actually represent the share of each cell’s area included in the shapefile (values from 0 to 1, i.e. 0 % to 100 %):
armenia_rastered <- rasterize(armenia, miroc_sample_C_heat, getCover=TRUE)
armenia_rastered[armenia_rastered==0] <- NA
plot(armenia_rastered)
plot(armenia, add=T)
We now use this newly created raster “armenia_rastered” as a mask for “miroc_sample_C_heat” and see that in the resulting layer “miroc_sample_C_heat_masked”, even cells that only marginally intersect with Armenia are included:
miroc_sample_C_heat_masked <- mask(miroc_sample_C_heat, armenia_rastered)
plot(miroc_sample_C_heat_masked)
plot(armenia, add=T)
To count the number of cells that are equal to one, we can just query the values contained in “miroc_sample_C_heat_masked”, i.e. the raster data itself:
sum(na.omit(values(miroc_sample_C_heat_masked)))## [1] 2Suppose we want to calculate the amount of heat events in January 2081. To avoid using the mask() function for each daily layer (which is computationally expensive), we can first calculate the sum of events in each cell for 31 days, and then perform the mask() function only once. The climate model predicts 83 heat events in Armenia for January 2081, i.e. 83 times a raster cell that intersects Armenia has a value above 10 degrees Celsius:
## select all layers of January 2081
january2081 <- miroc_stack_C[[1:31]]
## reclassify raster values
january2081[january2081 <= 10] <- 0
january2081[january2081 > 10] <- 1
## calculate sum
january2081_sum <- sum(january2081)
plot(january2081_sum)
plot(armenia, add=T)
## mask all raster cells overlapping Armenia
january2081_sum_masked <- mask(january2081_sum, armenia_rastered)
plot(january2081_sum_masked)
plot(armenia, add=T)
## count January heat events for whole Armenia
sum(na.omit(values(january2081_sum_masked)))## [1] 83If we want to select a different month in a different year, we can make use of the layer names of “miroc_stack_C” to identify the respective positions. The layer names are a vector of character values that all have exactly the same amount of characters, so we can use the substring() function to look for a certain pattern:
head(names(miroc_stack_C))## [1] "X2081.01.01" "X2081.01.02" "X2081.01.03" "X2081.01.04" "X2081.01.05"
## [6] "X2081.01.06"## identify the layer positions of a specific year:
which(substring(names(miroc_stack_C), 2, 5) == "2085")## [1] 1462 1463 1464 1465 1466 1467 1468 1469 1470 1471 1472 1473 1474 1475 1476
## [16] 1477 1478 1479 1480 1481 1482 1483 1484 1485 1486 1487 1488 1489 1490 1491
## [31] 1492 1493 1494 1495 1496 1497 1498 1499 1500 1501 1502 1503 1504 1505 1506
## [46] 1507 1508 1509 1510 1511 1512 1513 1514 1515 1516 1517 1518 1519 1520 1521
## [61] 1522 1523 1524 1525 1526 1527 1528 1529 1530 1531 1532 1533 1534 1535 1536
## [76] 1537 1538 1539 1540 1541 1542 1543 1544 1545 1546 1547 1548 1549 1550 1551
## [91] 1552 1553 1554 1555 1556 1557 1558 1559 1560 1561 1562 1563 1564 1565 1566
## [106] 1567 1568 1569 1570 1571 1572 1573 1574 1575 1576 1577 1578 1579 1580 1581
## [121] 1582 1583 1584 1585 1586 1587 1588 1589 1590 1591 1592 1593 1594 1595 1596
## [136] 1597 1598 1599 1600 1601 1602 1603 1604 1605 1606 1607 1608 1609 1610 1611
## [151] 1612 1613 1614 1615 1616 1617 1618 1619 1620 1621 1622 1623 1624 1625 1626
## [166] 1627 1628 1629 1630 1631 1632 1633 1634 1635 1636 1637 1638 1639 1640 1641
## [181] 1642 1643 1644 1645 1646 1647 1648 1649 1650 1651 1652 1653 1654 1655 1656
## [196] 1657 1658 1659 1660 1661 1662 1663 1664 1665 1666 1667 1668 1669 1670 1671
## [211] 1672 1673 1674 1675 1676 1677 1678 1679 1680 1681 1682 1683 1684 1685 1686
## [226] 1687 1688 1689 1690 1691 1692 1693 1694 1695 1696 1697 1698 1699 1700 1701
## [241] 1702 1703 1704 1705 1706 1707 1708 1709 1710 1711 1712 1713 1714 1715 1716
## [256] 1717 1718 1719 1720 1721 1722 1723 1724 1725 1726 1727 1728 1729 1730 1731
## [271] 1732 1733 1734 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746
## [286] 1747 1748 1749 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761
## [301] 1762 1763 1764 1765 1766 1767 1768 1769 1770 1771 1772 1773 1774 1775 1776
## [316] 1777 1778 1779 1780 1781 1782 1783 1784 1785 1786 1787 1788 1789 1790 1791
## [331] 1792 1793 1794 1795 1796 1797 1798 1799 1800 1801 1802 1803 1804 1805 1806
## [346] 1807 1808 1809 1810 1811 1812 1813 1814 1815 1816 1817 1818 1819 1820 1821
## [361] 1822 1823 1824 1825 1826## identify the layer positions of a specific year, and a specific month:
which( (substring(names(miroc_stack_C), 2, 5) == "2085") & (substring(names(miroc_stack_C), 7, 8) == "10") )## [1] 1735 1736 1737 1738 1739 1740 1741 1742 1743 1744 1745 1746 1747 1748 1749
## [16] 1750 1751 1752 1753 1754 1755 1756 1757 1758 1759 1760 1761 1762 1763 1764
## [31] 17653. Calculate extreme events using an own function
Let us put the steps outlined above into a function that returns the number of heat events in a given month and year, i.e. the number of cells that intersect Armenia and exceed a certain threshold during this period of time. This function needs four arguments - (i) the RasterStack object that should be processed (= “data”), (ii) the threshold to be exceeded (= “threshold”), (iii) the month of interest (= “month”) and (iv) the year of interest (= “year”). We call this function “CountCellsArmenia”:
CountCellsArmenia <- function(data, threshold, month, year) {
data_sel <- data[[which( (substring(names(data), 2, 5) == year) & (substring(names(data), 7, 8) == month) )]]
data_sel[data_sel <= threshold] <- 0
data_sel[data_sel > threshold] <- 1
data_sel_sum <- sum(data_sel)
data_sel_sum_masked <- mask(data_sel_sum, armenia_rastered)
return(sum(na.omit(values(data_sel_sum_masked))))
}Check if “CountCellsArmenia” yields the same result for January 2081 as before:
CountCellsArmenia(miroc_stack_C, 10, "01", 2081)## [1] 83Let’s iteratively run the function for all months of July from 2081 to 2090 and save the results in a vector. Let’s take a heat threshold of 40 degrees Celsius here, and plot the result as a barplot:
years <- c(2081:2090)
July_events <- data.frame(year=as.factor(years), heat_events=NA)
for (y in 1:length(years))
{ July_events$heat_events[y] <- CountCellsArmenia(miroc_stack_C, 40, "07", years[y]) }
July_events| year | heat_events |
|---|---|
| 2081 | 86 |
| 2082 | 28 |
| 2083 | 126 |
| 2084 | 223 |
| 2085 | 114 |
| 2086 | 35 |
| 2087 | 62 |
| 2088 | 22 |
| 2089 | 39 |
| 2090 | 133 |
ggplot(data=July_events, aes(x=year, y=heat_events)) +
geom_bar(stat="identity") +
labs(y="raster cells exceeding threshold", x="") +
theme(axis.text.x=element_text(angle=45, hjust=1)) +
ggtitle("Heat events in Armenia according to MIROC5 under RCP8.5. Threshold: 40 °C")
We do not see a trend here. However, note that in these future models, the exact values during a day or even during one year should not be over-interpreted. Usually, you take the average over 10 or 20 years to describe climate and weather patterns in the future. The average yearly amount of heat events in Armenia would here be 86.8 with a standard deviation of 63.4, and it would make sense to compare this number to other climate models, other RCP scenarios or other future decades:
mean(July_events$heat_events)## [1] 86.8sd(July_events$heat_events)## [1] 63.42239plotdata <- data.frame(month="July", mean=mean(July_events$heat_events), sd=sd(July_events$heat_events))
ggplot(data=plotdata, aes(x=month, y=mean)) +
geom_bar(stat="identity") +
geom_errorbar(aes(ymin=mean-sd, ymax=mean+sd)) +
labs(y="raster cells exceeding 40 °C", x="") +
ggtitle("Mean yearly heat events in Armenia according to MIROC5 under RCP8.5 from 2081 to 2090") 
Congratulations, you made it to the end of the second script! You can now solve the exercises 06 to 10.