2.3. Reading and writing spatio-temporal data¶

Python data types (e.g. floating points) can be printed, written to disk, or read from disk using functions that come with standard Python. Think of print to print a variable on the screen, or the file objects to read and write data to disk. As these functions do not work with PCRaster maps, PCRaster Python provides its own functions for reading and writing map data.

2.3.1. Reading and writing static spatial data¶

Open the script dynMod.py and save it under the name staticMap.py. In the initial, add the following lines of code:

dem = self.readmap('dem')
slopeOfDem = slope(dem)
self.report(slopeOfDem,'gradient')

Run the model script. The readmap is a function to open a map stored on disk, you do not need to give the suffix .map of the filename that is opened (dem.map). The report function writes a map (here, slopeOfDem) to disk under a name given as the second argument, which is of data type string. When used in the initial (as you did here), it stores a single map, here the filename will be gradient.map.

Check the content of gradient.map by displaying it.

2.3.2. Reading and writing temporal spatial data¶

A very similar approach can be followed for temporal map data. However now we need to add statements to the dynamic. Let’s import the time series of precipitation maps that you displayed in the visualisation part of the exercises to the model, convert it to mm/hours and write the result as a timeseries of maps.

Open staticMap.py created in the previous section and save it as dynamicMap.py. In the dynamic section:

Remove pass.
Add a statement that reads the precipitation from disk using self.readmap(`precip`).
Add a statement converting from m/day to mm/day.
Write the result to disk using self.report (syntax is the same as in the initial). As second argument of the report function, use pmm.

Also, change the number of timesteps that the model is run to 181.

Save and run the script. In the command shell, type dir (or ls on unix) to see what is stored. If everything is OK, you should see the filenames pmm00000.001, pmm00000.002, etc.

Question: Use aguila to animate the original precipitation timseries of maps that you read from disk, together with the same precipitation converted to mm/day (one Aguila command). Check the results. What command did you use?

aguila --timesteps=[1,181,1] precip pmm
aguila --timesteps=1,181,1 precip pmm
aguila -timesteps=[1,181,1] precip pmm
aguila --timesteps=[1,181,1] precip.tss pmm.tss

Now, try to make an animation of maps that shows where precipitation is above 0.01 m/day. Add a comparison operation to the dynamic that calculates a Boolean map with True where cells are above the threshold and False elsewhere. If needed use the PCRaster man pages. Store the output on harddisk using report. Run the model and animate the results with Aguila.

Question: At what elevation do we find particularly high values of precipitation?

There is no relation with elevation.
Large precipitation is mainly found in the valleys.
Precipitation is zero throughout the area.
Precipitation is high in the higher areas.

2.3.3. Reading timeseries¶

Timeseries files can be imported to a script using the timeinput.. function. As most PCRaster functions operate on maps, the timeinput.. function reads the value(s) from the time series for the current time step, and directly converts it to a PCRaster map that can be used in calculations. Let’s read from disk the precipitation timeseries file that you had a look at in the visualisation section of the exercises.

Open dynMod.py and save it as dynamicTimeSeries.py. Change the number of time steps to 181 and replace the pass in the dynamic with the line:

precipitation=timeinputscalar('precip.tss',1)

As the precip.tss does not contain information on the data type of the data in the file, you need to provide this with the timeinput command. Precipitation are scalar data, so we use timeinputscalar. In case of nominal data, for instance, you would use timeinputnominal. The first input argument is the timeseries that is read. The second argument is the area to which the values in the timeseries is assigned. Here we want to assign the values in the timeseries value to the whole map, so we provide a 1 as second argument. This should be read as a map completely filled with one values. In case the timeseries contains multiple columns, i.e. multiple locations, you need to enter a map with regions as second argument. See the PCRaster man page on timeinput.. for the details.

Add a report statement that writes precipitation to disk. Animate the results with Aguila.

Question: What is the content of the resulting timeseries?

The precipitation value from the .tss file is uniformly assigned to the whole map for each time step.
It contains a value one for each time step.
The precipitation is distributed according to the elevation.