In this use case example we will calculate the yearly revealing of dry land due to land uplift in Finland. This is done by constructing paleotopographies and calclulating the areal change between different time periods.
Even though the yearly speed of land uplift in Finland is relatively well studied, the amount of revealing land has been a bit problematic to calculate. With the current national data sets and high performance computing services of CSC, we are able to compute some estimates of the yearly land revealing.
The batymetric data of Finnish coastal areas is relatively inaccurate and unevenly distributed, this is why we approach the issue using the terrestrial elevation models. This means we will travel back in time and estimate the current state of land revealing by calculating the amount of land revealing in the past. We will do this by constructing paleotopographies and comparing them to each other. In the computing of paleotopographies we will make use of current terrestrial elevation models and land uplift model.
As the elevation models tend to contain errors and false information at the water areas, the sea areas will be masked to value 0 with the sea polygons of Topographic Database of Finland. Because the coastline of the Topographic Database doesn't always match with the digital elevation models, we will start our calculations from the year 200 bp when the coastline was located further inland than it is today. (before present). From there we will construct 10 paleotopographies between every 50 years finishing at year 700bp.
Current land uplift rates can be used untill the year 1890. From there on, we will add 1 mm to the yarly rates as suggested by Ekman (2001).
- 2m Elevation model by National Land Survey of Finland. Data available in Taito.
- 10m Elevation model by National Land Survey of Finland. Data available in Taito.
- UTM map sheet division by the National Land Survey of Finland. Data available in the File service of open data. The needed UTM_10 grid is stored for you in data folder.
- Sea areas of the Topographic Database by the National Land Survey of Finland. Data available in Taito
- Isostasy point data based on the NKG2016LU_lev land uplift data by the Nordic Geodetic Comission. The point data available in data folder is allready cropped to cover only the coastal areas of Finland. Read more about the data in here.
Scripts used : 10mDem_masker_resampler.py, dem10_batch
Because the 2m elevation model is not available everywhere in Finland, we need to use 10m elevation model at some areas. To make sure that we have the data when needed, we will prepare it in advance.
Run the script 10mDem_masker_resampler.py using the dem10_batch. At the first phase the script will go through the directory of 10m dem files, mask the sea areas to value 0 with the intersecting sea area polygons. The masking to 0 is essential because the water areas of the laser scanning data include a lot of errors. Negative values are also converted to 0. The masked file is saved if it fills the following criteria:
- We are only interested in dem files with minimum value less than 7.3. Higher dem files would have no changes in them. The value is based on the calculated maximum of land uplift in 700 years.
(700yr x 9,34mm) + ((1890yr-(2019yr-700yr)) x 1mm) = 7109mm
Where 9.34mm is the maximum value of yearly land uplift at the study area. Also the 1mm addition to land uplift values after year 1890 is taken into account (Ekman 2001). Safety marginal of 0.2m was added to the result.
- We are also only interested in dem files that have dry land in them. This is why the max value of the file has to be over 0m.
In the second phase the masked 10m dem files are resampled to 2m resolution. The 10m dem needs to be resampled to 2m resolution so we can make calculations between other layers. After resampling the files are renamed and saved.
Scripts used: dryland_calculator.py , calculator_batch
You can run the dryland_calculator using the calculator_batch file in Taito. The script is designed so that it runs as an array job. The jobs are separated with the first two letters of UTM10 grid ID. This way we will also limit the search to the coastal areas of Finland. The used map sheets and their ID.s are shown on the map.
The actual calculation of the land revealing is done in the UTM10 map division level, which means 6km x 6km grid squares that are named for example ”L4314A”. This is because the 2m dem is stored in files named and defined after this level of division.
The calculator goes through a file containing the ID ("LEHTITUNNU") and the geometry of UTM10 grid cells. Whenever a corresponding 2m dem file is found, it is masked and filtered just like 10m dem files previously. If a corresponding 2m dem file cannot be found, the calculator opens the already resampled and masked 10m dem file and uses that in the calculations instead of the 2m dem.
Before we can begin the calculations we also need a 2m resolution raster layer of the isostasy data. The isostasy raster is created on the fly with each iteration only to the area corresponding the dem file by using linear interpolation.
The first phase of the calculations between the dem layer and the isostasy layer is to construct a paleotopography for every time period. In this example we use 10 time periods between every 50 years from 200bp. to 700bp. The paleotopography can be calculated by multiplying the isostasy layer with the number of years and then extracting the isostasy layer form the present day dem file. Before extraction the isostasy values were converted to meters. The change in the speed of isostasy after 1890 was also taken into account. The ten paleotopography layers for each grid ID were calculated as following:
dem_year = dem - ((isostasylayer*year)+(year-(year-(2019-1890)))/1000)
In the next phase we want to calculate the difference in pixels that have values greater than 0 between consecutive paleotopographies. This way we can establish how many pixels were revealed during the 50 year time interval. Then we get the amount of revealed land by multiplying the number of changed pixels with their area in square kilometres (4m^2 / 1000000). The calculations are done for all the 10 time intervals as following:
change = ((0 < dem_year1).sum()-(0 < dem_year2).sum())*4/1000000
The 10 resulting values are saved in text file with the grid cell id and stored in folders based on the map sheet division. One result text file will look something like this:
0.084368,0.07966,0.075436,0.070012,0.071252,0.063976,0.059912,0.05324,0.046804,0.041728,K3222F
Script used: result_gatherer.py
Because all the results are stored in separate text files, we need to gather them together. Result gatherer will go through all the text files in the directory and gather the data into one table. This can be done for example at your local computer or in Taito through NoMachine Spyder interface.
After the results are gathered it is relatively easy to plot the results for example as bar plots:
You can also export an shapefile merging the utm10 grid with the data you just created. This way you can visualize maps for example in QGIS or ArcMap:
According to the results, the estimation of present yearly formation of dry land is approximately 8.9 km2. The yearly land revealing doesn't seem to be consistent and is increasing when we move towards present day. As expected, the revealing of new land is strongest at the western coast of Finland where also the speed of land uplift is greatest.
Ekman, M (2001) Calculation of Historical Shore Levels in Fennoscandia due to Postglacial Rebound. Small Publications in Historical Geophysics 8
When using the scripts or CSC.s computation services, please cite the oGIIR project: "We made use of geospatial data/instructions/computing resources provided by CSC and the Open Geospatial Information Infrastructure for Research (oGIIR, urn:nbn:fi:research-infras-2016072513) funded by the Academy of Finland."
Authored by Akseli Toikka and the Department of Geoinformatics and Cartography at FGI