LSStats

Screenshot

Download Model

Download the .scn, .sel and .lse files by clicking on the following link: File:LSStats.zip

Model Code Exploration

In the following sections we will examine all of the model files for this model. Note that instead of downloading the zip file above, you could just copy the text in the boxes below into a text editor and save it with the appropriate name (Section title). Opening the resulting .scn file in the SELES simulator would run this model.

BatchStats.scn

SELES Scenario
$baseGisDir$ = ..\..\CaseStudy\gisData\cell
$gisDir$ = ..\..\CaseStudy\v6_roads\output\cell
$first$ = TRUE
// load first layers
StudyArea = $baseGisDir$\StudyArea
cwd oStats
// Loop over all files in the directory that start with "ACPatchType1."
for($x$ = "..\$gisDir$\OldForest_1_*")
// Open a patch layer (ensure previous one, if any, is closed)
PatchLayer = ..\$gisDir$\OldForest_1_$x$
// Load model
Model Dimensions: PatchLayer
..\stats.sel
// Set up parameters
CellWidth = 100
HaPerCell = 1
NNType = 2         // For MPG
LandscapeId = $x$  // This sets the variable LandscapeId to the integer represented by $x$
  NumPTypes = 1
  Reset Output $first$     // Tell SELES to not reset output files (i.e. to append to them) except on first  iteration
  $first$ = FALSE
  // Turn off centroid stats (time consuming if there are lots of patches (np^2))
//   centroid.lse OFF
//   nn.lse OFF
  // Run simulation
  SimStart 100*365.25 1
  Close PatchLayer
end

centroid.lse

LSEVENT: CENTROID_LSE
DEFINITIONS
  GLOBAL CONSTANT: MaxPTypes
  GLOBAL VARIABLE: NumPTypes, MinPType, HaPerCell, CellWidth, LandscapeId, Replicate
  /* Assume these stats are defined */
  GLOBAL VARIABLE: A, Area[], NP, NumPatches[]
  /* Mean Connectivity between patches */
  GLOBAL VARIABLE: MeanCCE[MaxPTypes+1], MaxCCE[MaxPTypes+1], MeanCD[MaxPTypes+1]
  GLOBAL VARIABLE: tCCE, tMaxCCE, tCD
  // Patch list: type and area previously defined
  GLOBAL CONSTANT: NumPatchStats, rId, rType, rArea, rCentroid, rMaxCCE
  GLOBAL LIST{NumPatchStats} VARIABLE: patchList[]
  GLOBAL VARIABLE: patchVar[NumPatchStats]
  EVENT VARIABLE: id
  CLUSTER VARIABLE: rowTotal, colTotal, numActiveCells
//   MIN OUTPUT VARIABLE: LSStatsFile = "LSStatsC.txt"
  MIN OUTPUT VARIABLE: ClassStatsFile = "ClassStatsC.txt"
//   MIN OUTPUT VARIABLE: EdgeStatsFile = "EdgeStatsC.txt"
ENDDEF
INITIALSTATE
  tCCE = 0
  tMaxCCE = 0
  tCD = 0
  OVER INDEX SEQUENCE(0, NumPTypes-1)
     pType = Index + MinPType
     MeanCCE[pType] = 0
     MaxCCE[pType] = 0
     MeanCD[pType] = 0
  ENDFN
  /* All the stats area available in the patch list, we only need to */
  /* summarize and output results */
  INITIALSTATE = 1
ENDIS
RETURNTIME
//   RETURNTIME = IF Time EQ 0 THEN (4 * 365.25) ELSE 0
  RETURNTIME = IF Time EQ 0 THEN (4 * 365.25) ELSE 0
  /* Finalize stats */
  // Iterate over each pair of patches of the same type
  OVER INDEX SEQUENCE(0, NumPTypes-1)
     i = Index + MinPType
     pos1 = HEAD(patchList[i])
     WHILE (pos1)
        patchVar [=] GET(patchList[i], pos1)
        currType = patchVar[rType]
        id1 = patchVar[rId]
        patchArea1 = patchVar[rArea]
        centroid1 = patchVar[rCentroid]
        maxCCE1 = patchVar[rMaxCCE]
        nextPos = NEXT(patchList[i], pos1)
//         pos2 = HEAD(patchList[i])
        pos2 = nextPos
        WHILE (pos2)
           IF (currType EQ GET(patchList[i], pos2, rType))
              patchVar [=] GET(patchList[i], pos2)
              patchArea2 = patchVar[rArea]
              centroid2 = patchVar[rCentroid]
              maxCCE2 = patchVar[rMaxCCE]
             /* Determine the max patch CCE */
              d = DISTANCE(centroid1, centroid2)
              CCE = (patchArea1 * patchArea2) / d^2
              MeanCCE[currType] = MeanCCE[currType] + CCE
              maxCCE1 = MAX(CCE, maxCCE1)
              maxCCE2 = MAX(CCE, maxCCE2)
              MeanCD[currType] = MeanCD[currType] + d
              // Update maxCCE in patches
              SET(patchList[i], pos2, rMaxCCE, maxCCE2)
/*
              OUTPUT RECORD(EdgeStatsFile)
                 LandscapeId: LandscapeId
                 Replicate: Replicate
                 type: currType
                 patchId1: id1
                 patchId2: patchVar[rId]
                 d: d
              ENDFN
*/
           ENDFN
           pos2 = NEXT(patchList[i], pos2)
        ENDFN
       /* At this point, the patch CCE for p1 will be known */
        SET(patchList[i], pos1, rMaxCCE, maxCCE1)
        MaxCCE[currType] = MaxCCE[currType] + maxCCE1
        pos1 = nextPos
     ENDFN
  ENDFN
  OVER INDEX SEQUENCE(0, NumPTypes-1)
     currType = Index + MinPType
     tCCE = tCCE + MeanCCE[currType]
     tMaxCCE = tMaxCCE + MaxCCE[currType]
     MeanCCE[currType] = MeanCCE[currType] / (NumPatches[currType] * (NumPatches[currType] - 1) * 0.5)
     MaxCCE[currType] = MaxCCE[currType] / (NumPatches[currType] - 1)
     tCD = tCD + MeanCD[currType]
     MeanCD[currType] = MeanCD[currType] / (NumPatches[currType] * (NumPatches[currType] - 1) * 0.5)
     OUTPUT RECORD(ClassStatsFile)
        DECISION Area[currType] > 0
        LandscapeId: LandscapeId
        Replicate: Replicate
        pType: currType
        CCE: MeanCCE[currType]
        MaxCCE: MaxCCE[currType]
        CD: MeanCD[currType]
     ENDFN
  ENDFN
  tCCE = IF (NP <= 1) THEN 0 ELSE tCCE / (NP * (NP - 1) * 0.5)
  tMaxCCE = tMaxCCE / NP
  tCD = IF (NP <= 1) THEN 0 ELSE tCD / (NP * (NP - 1) * 0.5)
/*
  OUTPUT RECORD(LSStatsFile)
     LandscapeId: LandscapeId
     Replicate: Replicate
     tCCE: tCCE
     tMaxCCE: tMaxCCE
     tCD: tCD
  ENDFN
*/
ENDRT
EVENTOPENINGS = 0

identifyPatches.lse

/* Combination of stats from fragstats and apack */
/* Note: A particular problem should decide on stats required */
/* and adapt a custom stats event, which may include some of the */
/* stats herein, and may include new or modified stats. */
/* Apack stats are denoted using lower case identifiers. */
LSEVENT: STATS_LSE
DEFINITIONS
  LAYER: PatchLayer, StudyArea, PatchId, CoreAreaId
  LAYER: EdgeInterior
    FEATURE 1: Edge
    FEATURE 2: Interior   
  GLOBAL CONSTANT: MaxPTypes
  CONSTANT: pi = 3.141593
  GLOBAL VARIABLE: NumPTypes, MinPType, HaPerCell, CellWidth, LandscapeId
     // Area indices
  GLOBAL VARIABLE: Area[]
  // Patch indices
  GLOBAL VARIABLE: NumPatches[]
  // Edge Indices
  GLOBAL VARIABLE: TotalEdge[], TotalEdgeArea[]
  // Core Area Indices
  GLOBAL VARIABLE: TCA[], NumCA[]
  // Adjacency matrix
  GLOBAL VARIABLE: Eik[]
  // pContagion
  GLOBAL VARIABLE: NumNeighbs[], NumSameNeighbs[]
  // Patch list
  GLOBAL CONSTANT: NumPatchStats, rId, rType, rArea, rCoreArea, rPerim, rCentroid, rMaxCCE, rNNDist, MaxDist
  GLOBAL LIST{NumPatchStats} VARIABLE: patchList[], coreAreaList[]
  //   GLOBAL LIST{NumPatchStats} VARIABLE: patchList, coreAreaList
  GLOBAL VARIABLE: patchVar[NumPatchStats]
  GLOBAL VARIABLE: pId, coreId
  EVENT VARIABLE: isPatchEvent, currType
  CLUSTER VARIABLE: numActiveCells, currPatchSize, currPatchCoreSize, currPatchPerim
  CLUSTER VARIABLE: rowTotal, colTotal
ENDDEF
INITIALSTATE
  EdgeInterior = 0
  pId = 0
  coreId = 0
  /* Three types of instances: the first NumPType instances identify and computes first order stats for patches */
  /*  (a separate instance is used per type so that the list is ordered by type) */
  /* The second identifies  and computes first order stats for core areas */
  /* The third computes second order stats and then summarizes and outputs */
  INITIALSTATE = 1 + NumPTypes
  isPatchEvent = EventId <= NumPTypes
  currType = EventId + MinPType - 1   // This is ok for patches, but don't use for core areas
ENDIS
RETURNTIME = IF Time EQ 0 THEN EventId + 365.25 ELSE 0
PROBINIT
  PROBINIT = IF (isPatchEvent) THEN (PatchLayer EQ currType) AND (StudyArea > 0)
             ELSE (EdgeInterior EQ Interior)
  numActiveCells = 0
  currPatchSize = 0
  currPatchCoreSize = 0
  currPatchPerim = 0
  rowTotal = 0
  colTotal = 0
  // For the patches ...
  IF (isPatchEvent)
     Area[PatchLayer] = Area[PatchLayer] + 1
     IF (PatchId EQ 0) // not visited...
        pId = pId + 1
        NumPatches[PatchLayer] = NumPatches[PatchLayer] + 1
     ENDFN
  // For the core areas ...
  ELSE
     coreId = coreId + (CoreAreaId EQ 0)
     NumCA[PatchLayer] = NumCA[PatchLayer] + (CoreAreaId EQ 0)
  ENDFN
ENDPI
TRANSITIONS
  TRANSITIONS = IF isPatchEvent THEN (PatchId EQ 0) ELSE (CoreAreaId EQ 0)
  // For the patches ...
  IF isPatchEvent
     numSimilarNeighbours = 0
     numDifferentNeighbours = 0
     OVER REGION CENTRED(1, 1)
        DECISION (0 <= PatchLayer < (MinPType + NumPTypes)) AND (StudyArea > 0)
        similarNeighb = (PatchLayer EQ currType)
        numSimilarNeighbours = numSimilarNeighbours + similarNeighb
        numDifferentNeighbours = numDifferentNeighbours + (similarNeighb EQ FALSE)
        NumNeighbs[PatchLayer] = NumNeighbs[PatchLayer] + 1
        NumSameNeighbs[PatchLayer] = NumSameNeighbs[PatchLayer] + similarNeighb 
        Eik[currType, PatchLayer] = Eik[currType, PatchLayer] + (similarNeighb EQ FALSE)
     ENDFN
     currPatchPerim = currPatchPerim + 4 - numSimilarNeighbours
     TotalEdge[PatchLayer] = TotalEdge[PatchLayer] + numDifferentNeighbours
     // Add in the diagonal neighbours
     OVER REGION CENTRED(1.1, 1.5)
        DECISION (0 <= PatchLayer < (MinPType + NumPTypes))
        similarNeighb = (PatchLayer EQ currType)
        numSimilarNeighbours = numSimilarNeighbours + similarNeighb 
        numDifferentNeighbours = numDifferentNeighbours + (similarNeighb EQ FALSE)
     ENDFN
     EdgeInterior = IF (numDifferentNeighbours > 0) THEN Edge ELSE Interior
     TCA[PatchLayer] = TCA[PatchLayer] + (EdgeInterior EQ Interior)
     TotalEdgeArea[PatchLayer] = TotalEdgeArea[PatchLayer] + (EdgeInterior EQ Edge)
     PatchId = pId
  ELSE
     CoreAreaId = coreId
  ENDFN
  currPatchSize = currPatchSize + 1
  currPatchCoreSize = currPatchCoreSize + (EdgeInterior EQ Interior)
  numActiveCells = numActiveCells + 1
  rowTotal = rowTotal + ROW(Location)
  colTotal = colTotal + COL(Location)
ENDTR
SPREADLOCATION
  REGION CENTRED(1, 1.5)
     DECISION IF isPatchEvent THEN (PatchId EQ 0) AND (PatchLayer EQ currType) AND (StudyArea > 0)
              ELSE (CoreAreaId EQ 0) AND (EdgeInterior EQ Interior) //AND (PatchLayer EQ currType)
ENDSL
SPREADTIMESTEP
  SPREADTIMESTEP = -1
  // If numActiveCells becomes 0, then this patch is done
  numActiveCells = numActiveCells - 1
  IF (numActiveCells EQ 0)
     // Add the patch to the list
     patchVar[rArea] = currPatchSize
     patchVar[rCoreArea] = currPatchCoreSize
     patchVar[rPerim] = currPatchPerim
     patchVar[rCentroid] = LOCATION(FLOOR(rowTotal/currPatchSize), FLOOR(colTotal/currPatchSize))
     patchVar[rMaxCCE] = 0
     patchVar[rNNDist] = -1
     IF isPatchEvent
        patchVar[rType] = currType
        patchVar[rId] = pId
        INSERT TAIL(patchList[currType], patchVar)
     ELSE
        patchVar[rType] = PatchLayer
        patchVar[rId] = coreId
        INSERT TAIL(coreAreaList[PatchLayer], patchVar)
     ENDFN
  ENDFN
ENDST

stats.sel

Seles Model
Time Units: Step kiloStep 1000 100000
Landscape Events:
ConvexHull.lse
Spatial Constants:
StudyArea
Spatial Variables:
ConvexHullPoints[1] <= 0
ConvexHull[1] <= 0
Global Variables:
TRLoc = 0
Output Frequency: 1

ConvexHull.lse

// Sweep-line approach to identifying convex hull using a worker agent
LSAGENT: ConvexHull
DEFINITIONS
 LAYER: StudyArea, ConvexHull, ConvexHullPoints
 GLOBAL VARIABLE: TRLoc
 LOCAL VARIABLE: EndingPivotLoc
 AGENT VARIABLE: PivotLocation, EdgeLocation, EdgeRow, EdgeCol
ENDDEF
INITIALSTATE
  INITIALSTATE = 1
// Find the top-right-most point
  TRLoc = 0
  OVER REGION WHOLE MAP
     DECISION StudyArea > 0
     TRLoc = MAX(TRLoc, Location)
  ENDFN
ConvexHull = 1 // Initialize convex hull layer to all ones (will erase areas outside hull)
ENDIS
NUMAGENTS = 1 // Only need a single worker
// Start on the first pivot vertex on convex hull
AGENTLOCATION
  REGION LOCATION(TRLoc)
ENDAL
PROBINIT
 PROBINIT = 1
// Set up first pivot vertex, and first raster boundary location
 PivotLocation = Location
 EdgeRow = NUMROWS-1
 EdgeCol = COL(Location)
 EdgeLocation = LOCATION(EdgeRow, EdgeCol)
 EndingPivotLoc = -1
ENDPI
TRANSITIONS
 // Continue until we reach the ending pivot location (which is actually second time we hit the second pivot           vertex,
 // since the first pivot vertex gets only partially processed on first pass)
 TRANSITIONS = (PivotLocation NEQ EndingPivotLoc)
ConvexHullPoints = 1
// Walk line from boundary to vertex and erase until a new pivot vertex is hit (then move there)
 oldPivot = PivotLocation
 hitPivot = FALSE
 OVER REGION VECTOR(EdgeLocation, PivotLocation)
    DECISION !hitPivot
 IF (StudyArea > 0) // hit new vertex
       hitPivot = TRUE
       IF (EndingPivotLoc EQ -1) AND (Location NEQ PivotLocation)
          EndingPivotLoc = PivotLocation
       ENDFN
       PivotLocation = Location
    ELSE
       ConvexHull = 0 // erase
    ENDFN
 ENDFN
// Move boundary location clockwise around outside edge of raster
 IF (EdgeRow EQ NUMROWS-1) AND (EdgeCol < NUMCOLS-1)
    EdgeCol = EdgeCol + 1
 ELSE IF (EdgeCol EQ NUMCOLS-1) AND (EdgeRow > 0)
    EdgeRow = EdgeRow - 1
 ELSE IF (EdgeRow EQ 0) AND (EdgeCol > 0)
    EdgeCol = EdgeCol - 1
 ELSE
    EdgeRow = EdgeRow + 1
 ENDFN
 EdgeLocation = LOCATION(EdgeRow, EdgeCol)
ENDTR
POPULATIONTIME = 1
// Move to current/new pivot location
MOVELOCATION
 REGION LOCATION(PivotLocation)
ENDML

Suggested Experiments

To explore this cellular automata model further, try the following: