ForEST Demo 2

This ForEST demo will walk through a simple raster example using ASCII Grid files. Feel free to check out ForEST Demo 1 for additional tricks using map algebra in ForEST.


In [19]:
# Let's pull ForEST in using import
from forest import *

In [20]:
# A choice of Engines again (Multi-Processing, Tile, or the Pass (do nothing) Engine)
Config.engine=mp_engine
Config.engine=tile_engine
Config.engine=pass_engine

In [21]:
# This function will create some raster files in the ASCII Grid format and store them in the unittests folder.

def makefiles(nfiles,nrows,ncols):
   for f_i in range(nfiles):
        f = open("unittests/tmp_raster"+str(f_i)+".asc","w")
        f.write("ncols "+str(ncols)+"\n")
        f.write("nrows "+str(nrows)+"\n")
        f.write("xllcorner 0.0\n")
        f.write("yllcorner 0.0\n")
        f.write("cellsize 1.0\n")
        f.write("NODATA_value -999\n")
        
        for i in range(nrows):
            for j in range(ncols):
                f.write(str(i+j+f_i)+" ")
            f.write("\n")

        f.close()

# Call the function to make 3 raster files of dimension 6x4 (nice and small for demonstration purposes)
makefiles(3,6,4)

In [22]:
# Now we will take the file names from each raster file

r1_name = "unittests/tmp_raster0.asc"
r2_name = "unittests/tmp_raster1.asc"
r3_name = "unittests/tmp_raster2.asc"
ro_name = "unittests/tmp_rastero.asc"

In [23]:
# Now the ForEST magic. # This is a line of code in the ForEST langauge.

run_primitive(AGLoad.file(r1_name) == AGLoad.file(r2_name) == AGLoad.file(r3_name) == RasterAdd == RasterSub == AGStore.file(ro_name))

# Definitions:
# run_primitive: This will execute the Primitives in the Pattern
# AGLoad: This will read in (load) an Ascii Grid (AG) file
# RasterAdd: This will add two raster datasets
# RasterSub: This will subtract two raster datasets
# AGStore: This will write out (store) an Ascii Grid (AG) file
# ==: This is called a 'Sequence' it means keep going or "next instruction"

# Pattern explanation
# The first three primitives (AGLoad) will push the 3 raster datasets onto the data stack
# RasterAdd: pop off the top two (r2, r3), add them together, and push the output to the stack
# RasterSub: pop off the top two (r1, RasterAdd output), subtract them, and push the output to the stack
# AGStore:   pop off the remaining raster and write it to unittests/tmp_rastero.asc

# This style of programming is called postfix notation (or Reverse Polish Notation)
# This was all the rage when pocket protectors and hand-held Hewlett-Packard calculators were hip
# It is also a very powerful way to program Stack Machines, which is the computational model for ForEST
# Postfix notation translates A + B to A B +
# Notice the AGLoads (e.g., A, B) happen before the operations (e.g., RasterAdd, RasterSub) in our example

# This line of ForEST is the same as the following lines of pseudo code
# r1 = AsciiGridRead(r1_name)
# r2 = AsciiGridRead(r2_name)
# r3 = AsciiGridRead(r3_name)

# ri = RasterAdd(r2, r3)
# ro = RasterSub(ro, r1)
# AsciiGridWrite(ro, ro_name)

# While it looks a little funny at first
# It is a very simple and very powerful system to enable scalable spatial-temporal computation.
# If you want to learn more, then look at our paper in the Proceedings of GeoComputation 2019.
# https://auckland.figshare.com/articles/Space-Time_is_the_Key_For_Expressing_Spatial-Temporal_Computing/9870416


load file unittests/tmp_raster0.asc
load file unittests/tmp_raster1.asc
AsciiGridLoad == (Sequence) AsciiGridLoad
Running AsciiGridLoad
Open unittests/tmp_raster1.asc
4 6 0.0 0.0 1.0 -999.0
4.0 4.0
bobdata [[1. 2. 3. 4.]
 [2. 3. 4. 5.]
 [3. 4. 5. 6.]
 [4. 5. 6. 7.]
 [5. 6. 7. 8.]
 [6. 7. 8. 9.]]
PUSH Bob (0.000000,0.000000) [4.000000,6.000000]
load file unittests/tmp_raster2.asc
AsciiGridLoad == (Sequence) AsciiGridLoad
Running AsciiGridLoad
Open unittests/tmp_raster2.asc
4 6 0.0 0.0 1.0 -999.0
4.0 4.0
bobdata [[ 2.  3.  4.  5.]
 [ 3.  4.  5.  6.]
 [ 4.  5.  6.  7.]
 [ 5.  6.  7.  8.]
 [ 6.  7.  8.  9.]
 [ 7.  8.  9. 10.]]
PUSH Bob (0.000000,0.000000) [4.000000,6.000000]
AsciiGridLoad == (Sequence) RasterAdd
Running AsciiGridLoad
Open unittests/tmp_raster2.asc
4 6 0.0 0.0 1.0 -999.0
4.0 4.0
bobdata [[ 2.  3.  4.  5.]
 [ 3.  4.  5.  6.]
 [ 4.  5.  6.  7.]
 [ 5.  6.  7.  8.]
 [ 6.  7.  8.  9.]
 [ 7.  8.  9. 10.]]
PUSH Bob (0.000000,0.000000) [4.000000,6.000000]
RasterAdd == (Sequence) RasterSub
Running RasterAdd
POP
POP
PUSH Bob (0.000000,0.000000) [4.000000,6.000000]
store file unittests/tmp_rastero.asc
RasterSub == (Sequence) AsciiGridStore
Running RasterSub
POP
POP
PUSH Bob (0.000000,0.000000) [4.000000,6.000000]
Running AsciiGridStore
POP
Open for writing  unittests/tmp_rastero.asc
bobdata [[ 3.  4.  5.  6.]
 [ 4.  5.  6.  7.]
 [ 5.  6.  7.  8.]
 [ 6.  7.  8.  9.]
 [ 7.  8.  9. 10.]
 [ 8.  9. 10. 11.]]