[GRASS-SVN] r65541 - in grass-addons/grass7/raster/r.futures: . r.futures.calib r.futures.devpressure r.futures.pga
svn_grass at osgeo.org
svn_grass at osgeo.org
Mon Jul 6 11:56:57 PDT 2015
Author: annakrat
Date: 2015-07-06 11:56:57 -0700 (Mon, 06 Jul 2015)
New Revision: 65541
Added:
grass-addons/grass7/raster/r.futures/FUTURES_inputs_diagram.png
grass-addons/grass7/raster/r.futures/Makefile
grass-addons/grass7/raster/r.futures/r.futures.calib/
grass-addons/grass7/raster/r.futures/r.futures.calib/Makefile
grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html
grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.py
grass-addons/grass7/raster/r.futures/r.futures.devpressure/
grass-addons/grass7/raster/r.futures/r.futures.devpressure/Makefile
grass-addons/grass7/raster/r.futures/r.futures.devpressure/devpressure_example.png
grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html
grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.py
grass-addons/grass7/raster/r.futures/r.futures.html
grass-addons/grass7/raster/r.futures/r.futures.pga/
grass-addons/grass7/raster/r.futures/r.futures.pga/Makefile
grass-addons/grass7/raster/r.futures/r.futures.pga/incentive.png
grass-addons/grass7/raster/r.futures/r.futures.pga/main.cpp
grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html
grass-addons/grass7/raster/r.futures/r.futures.pga/r_futures.png
grass-addons/grass7/raster/r.futures/r.futures.pga/r_futures_detail.png
Removed:
grass-addons/grass7/raster/r.futures/Makefile
grass-addons/grass7/raster/r.futures/main.cpp
grass-addons/grass7/raster/r.futures/r.futures.html
grass-addons/grass7/raster/r.futures/r_futures.png
grass-addons/grass7/raster/r.futures/r_futures_detail.png
Log:
r.futures: updated version of FUTURES model with new tools for data preparation
Added: grass-addons/grass7/raster/r.futures/FUTURES_inputs_diagram.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.futures/FUTURES_inputs_diagram.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Deleted: grass-addons/grass7/raster/r.futures/Makefile
===================================================================
--- grass-addons/grass7/raster/r.futures/Makefile 2015-07-06 18:23:20 UTC (rev 65540)
+++ grass-addons/grass7/raster/r.futures/Makefile 2015-07-06 18:56:57 UTC (rev 65541)
@@ -1,16 +0,0 @@
-MODULE_TOPDIR = ../..
-
-PGM = r.futures
-
-LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
-DEPENDENCIES = $(RASTERDEP) $(GISDEP)
-
-include $(MODULE_TOPDIR)/include/Make/Module.make
-
-EXTRA_CFLAGS = -Wall
-
-LINK = $(CXX)
-
-ifneq ($(strip $(CXX)),)
-default: cmd
-endif
Added: grass-addons/grass7/raster/r.futures/Makefile
===================================================================
--- grass-addons/grass7/raster/r.futures/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.futures/Makefile 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,14 @@
+MODULE_TOPDIR =../..
+
+PGM = r.futures
+
+SUBDIRS = \
+ r.futures.calib \
+ r.futures.devpressure \
+ r.futures.pga
+
+include $(MODULE_TOPDIR)/include/Make/Dir.make
+
+default: parsubdirs htmldir
+
+install: installsubdirs
Deleted: grass-addons/grass7/raster/r.futures/main.cpp
===================================================================
--- grass-addons/grass7/raster/r.futures/main.cpp 2015-07-06 18:23:20 UTC (rev 65540)
+++ grass-addons/grass7/raster/r.futures/main.cpp 2015-07-06 18:56:57 UTC (rev 65541)
@@ -1,2754 +0,0 @@
-
-/****************************************************************************
- *
- * MODULE: r.futures
- * AUTHOR(S): Ross K. Meentemeyer
- * Wenwu Tang
- * Monica A. Dorning
- * John B. Vogler
- * Nik J. Cunniffe
- * Douglas A. Shoemaker
- * Jennifer A. Koch
- * Vaclav Petras <wenzeslaus gmail com>
- * (See the manual page for details and references.)
- *
- * PURPOSE: Simulation of urban-rural landscape structure (FUTURES model)
- *
- * COPYRIGHT: (C) 2013-2014 by Meentemeyer et al.
- *
- * This program is free software: you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation, either version 3
- * of the License, or (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public
- * License along with this program. If not, see
- * <http://www.gnu.org/licenses/> or read the file COPYING that
- * comes with GRASS GIS for details.
- *
- *****************************************************************************/
-
-/**
- \file main.c
-
- The main file containing both the model code and the data handing part.
-
- A lot of parts of the code are dead code. May be useful as an alternative
- algorithms. Some of the dead code is commented out and might need to
- be revived if other parts will show as useful.
-
- The language of the code is subject to change. The goal is to use either
- C or C++, not both mixed as it is now.
-
- Major refactoring of the code is expected.
-*/
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <string.h>
-#include <math.h>
-#include <sys/time.h>
-#include <iostream>
-#include <fstream>
-#include <vector>
-
-extern "C"
-{
-#include <grass/gis.h>
-#include <grass/raster.h>
-#include <grass/glocale.h>
-}
-
-//#include "distance.h"
-
-/** used to flag cells that are off the lattice */
-#define _CELL_OUT_OF_RANGE -2
-
-/** used to flag cells that are not in this county */
-#define _CELL_OUT_OF_COUNTY -1
-
-/** used to flag cells that are valid */
-#define _CELL_VALID 1
-
-/** need to dynamically allocate this as putting all on stack will crash most compilers */
-#define _N_MAX_DYNAMIC_BUFF_LEN (1024*1024)
-#define _N_MAX_FILENAME_LEN 1024
-
-/** Max length of input line */
-#define N_MAXREADINLEN 8192
-
-/** for working out where to start reading data */
-#define _GIS_HEADER_LENGTH 6
-
-/** strictly-speaking, should probably parse these from GIS files */
-#define _GIS_NO_DATA_STRING "-9999"
-
-/** Value used to represent NULL (NoData) in some parts of the model */
-#define _GIS_NO_DATA_INT -9999
-
-/** use this for tDeveloped if cell still undeveloped */
-#define _N_NOT_YET_DEVELOPED -1
-
-/* algorithm to use */
-
-/** deterministic model */
-#define _N_ALGORITHM_DETERMINISTIC 1
-
-/** stochastic model with downweighting of logit probabilities and of devPressure over time */
-#define _N_ALGORITHM_STOCHASTIC_I 2
-#define _N_ALGORITHM_STOCHASTIC_II 3
-
-#define _MAX_RAND_FIND_SEED_FACTOR 25
-#define maxNumAddVariables 6
-
-/** maximal number of counties allowed */
-#define MAXNUM_COUNTY 50
-#define MAX_YEARS 100
-/// maximum array size for undev cells (maximum: 1840198 for a county within 16 counties)
-#define MAX_UNDEV_SIZE 2000000
-using namespace std;
-
-
-/* Wenwu Tang */
-/// use absolute paths
-char dirName[100];
-
-/// string for temporarily storage
-char tempStr[100];
-
-
-
-typedef struct
-{
-
- /** see #define's starting _CELL_ above */
- int nCellType;
-
- /** attraction to employment base; static */
- double employAttraction;
-
- /** distance to interchange; static */
- double interchangeDistance;
-
- /** road density; static */
- double roadDensity;
-
- /** x position on the lattice; static */
- int thisX;
-
- /** y position on the lattice; static */
- int thisY;
-
- /** whether this site is still undeveloped; varies. Note if bUndeveloped = 0 the values of employAttraction etc. are not to be trusted */
- int bUndeveloped;
- // TODO: is the following int or double (defined as double but usually casting to int in assignment, was as int in some print)
-
- /** development pressure; varies */
- double devPressure;
-
- /** stores whether this site has been considered for development, to handle "protection" */
- int bUntouched;
-
- /** timestep on which developed (0 for developed at start, _N_NOT_YET_DEVELOPED for not developed yet ) */
- int tDeveloped;
-
- /** multiplicative factor on the probabilities */
- double consWeight;
-
- /** additional variables */
- double additionVariable[maxNumAddVariables];
- int index_region;
- float devProba;
-} t_Cell;
-
-typedef struct
-{
-
- /** id of this cell */
- int cellID;
-
- /** value of the logit */
- double logitVal;
-
- /** whether or not cell previously considered...need to use this in the sort */
- int bUntouched;
-
- /** to support random pick based on their logitVal */
- double cumulProb;
-} t_Undev;
-
-typedef struct
-{
-
- /** array of cells (see posFromXY for how they are indexed) */
- t_Cell *asCells;
-
- /** number of columns in the grid */
- int maxX;
-
- /** number of rows in the grid */
- int maxY;
-
- /** total number of cells */
- int totalCells;
-
- /** array of information on undeveloped cells */
- t_Undev *asUndev;
-
- /** number of cells which have not yet converted */
- int undevSites;
-
- /** posterior sample from parcel size distribution */
- int *aParcelSizes;
-
- /** number in that sample */
- int parcelSizes;
- std::vector < std::vector < t_Undev > >asUndevs; //WT
- int num_undevSites[MAXNUM_COUNTY]; //WT
-} t_Landscape;
-
-typedef struct
-{
-
- /** size of the grids */
- int xSize;
-
- /** size of the grids */
- int ySize;
-
- /** file containing information on how many cells to transition and when */
- char *controlFile;
-
- /** files containing the information to read in */
- char *employAttractionFile;
- char *interchangeDistanceFile;
- char *roadDensityFile;
- char *undevelopedFile;
- char *devPressureFile;
- char *consWeightFile;
- char *probLookupFile;
- int nProbLookup;
- double *adProbLookup;
- // parameters provided in other ways
- // /* parameters that go into regression formula */
- // double dIntercept; /* 0.038884 */
- // double dEmployAttraction; /* -0.0000091946 */
- // double dInterchangeDistance; /* 0.000042021 */
- // double dRoadDensity; /* 0.00065813 */
- // double dDevPressure; /* -0.026190 */
- /* size of square used to recalculate development pressure */
-
- int nDevNeighbourhood; /** 8 (i.e. 8 in each direction, leading to 17 x 17 = 289 in total) */
-
- /** used in writing rasters */
- char *dumpFile;
- char *outputSeries;
-
- /** 1 deterministic, 2 old stochastic, 3 new stochastic */
- int nAlgorithm;
-
- /** use this to downweight probabilities */
- double dProbWeight;
-
- /** and this to downweight old dev pressure */
- double dDevPersistence;
-
- /** file containing parcel size information */
- char *parcelSizeFile;
- double discountFactor; ///< for calibrate patch size
-
- /** give up spiralling around when examined this number of times too many cells */
- double giveUpRatio;
-
- /** these parameters only relevant for new stochastic algorithm */
- int sortProbs;
- double patchFactor;
- double patchMean;
- double patchRange;
- /// 4 or 8 neighbors only
- int numNeighbors;
- /// 1: uniform distribution 2: based on dev. proba.
- int seedSearch;
- int numAddVariables;
-
- /** parameters for additional variables */
- double addParameters[maxNumAddVariables][MAXNUM_COUNTY];
- char *addVariableFile[maxNumAddVariables];
- /// 1: #occurrence; 2: gravity (need to define alpha and scaling factor); 3: kernel, need to define alpha and scaling factor
- ///
- /// for 2: formula-> scalingFactor/power(dist,alpha)
- /// for 3: formula-> scalingFactor*exp(-2*dist/alpha)
- int devPressureApproach;
- /// scale distance-based force
- double scalingFactor;
- ///< constraint on distance
- double alpha;
-
- /** parameters that go into regression formula */
- double dIntercepts[MAXNUM_COUNTY];
- double dV1[MAXNUM_COUNTY];
- double dV2[MAXNUM_COUNTY];
- double dV3[MAXNUM_COUNTY];
- double dV4[MAXNUM_COUNTY];
- int num_Regions;
-
- /** index file to run multiple regions */
- char *indexFile;
-
- /// control files for all regions
- /// development demand for multiple regions
- char *controlFileAll;
- int devDemand[MAX_YEARS];
- int devDemands[MAXNUM_COUNTY][MAX_YEARS];
- /// number of simulation steps
- int nSteps;
-} t_Params;
-
-
-//From WT: move here to be global variables
-t_Params sParams;
-t_Landscape sLandscape;
-
-int getUnDevIndex(t_Landscape * pLandscape);
-void Shuffle4(int *mat);
-double getDistance1(double x1, double y1, double x2, double y2);
-void readDevPotParams(t_Params * pParams, char *fn);
-void readIndexData(t_Landscape * pLandscape, t_Params * pParams);
-void findAndSortProbsAll(t_Landscape * pLandscape, t_Params * pParams,
- int step);
-void updateMap1(t_Landscape * pLandscape, t_Params * pParams, int step,
- int regionID);
-int getUnDevIndex1(t_Landscape * pLandscape, int regionID);
-
-
-void readDevPotParams(t_Params * pParams, char *fn)
-{
- char str[200];
- int id;
- double di, d1, d2, d3, d4, val;
-
- // TODO: d5 is in file but unused
- int i, j;
- ifstream f;
-
- f.open(fn);
- if (!f.is_open())
- G_fatal_error(_("Cannot open development potential parameters file <%s>"),
- fn);
-
- f.getline(str, 200);
- for (i = 0; i < pParams->num_Regions; i++) {
- // TODO: read by lines to count the variables
- f >> id >> di >> d1 >> d2 >> d3 >> d4;
- cout << id << "\t" << di << "\t" << d1 << "\t" << d2 << "\t" << d3 <<
- "\t" << d4;
-
- pParams->dIntercepts[i] = di;
- pParams->dV1[i] = d1;
- pParams->dV2[i] = d2;
- pParams->dV3[i] = d3;
- pParams->dV4[i] = d4;
- for (j = 0; j < pParams->numAddVariables; j++) {
- f >> val;
- cout << "\t" << val;
- pParams->addParameters[j][i] = val;
- }
- cout << endl;
- }
- f.close();
-}
-
-/**
- Generate uniform number on [0,1).
-
- Encapsulated to allow easy replacement if necessary.
-*/
-double uniformRandom()
-{
- int nRet;
-
- nRet = RAND_MAX;
- while (nRet == RAND_MAX) { /* make sure never get 1.0 */
- nRet = rand();
- }
- return ((double)nRet / (double)(RAND_MAX));
-}
-
-/**
- Seed random number generator.
-
- Encapsulated to allow easy replacement if necessary.
-*/
-void seedRandom(struct timeval ttime)
-{
- srand((ttime.tv_sec * 100) + (ttime.tv_usec / 100));
- //srand((unsigned int)time(NULL));
-}
-
-
-/**
- Work out x and y position on the grid from index into siteMapping/landscape array.
-*/
-void xyFromPos(int nPos, int *pnX, int *pnY, t_Landscape * pLandscape)
-{
- /* integer division just gives whole number part */
- *pnX = nPos % pLandscape->maxX;
- *pnY = nPos / pLandscape->maxX;
-}
-
-/**
- Work out index into asCells array from location.
-*/
-int posFromXY(int nX, int nY, t_Landscape * pLandscape)
-{
- int nPos;
-
- /*
- 0 1 2 .... nX-1
- nX nX+1 nX+2 .... 2*nX-1
- . .
- . .
- . .
- (nY-1)*nX ... nX*nY-1
- */
- nPos = _CELL_OUT_OF_RANGE;
- if (nX >= 0 && nX < pLandscape->maxX && nY >= 0 && nY < pLandscape->maxY) {
- nPos = nX * pLandscape->maxY + nY;
- }
- return nPos;
-}
-
-/**
- Allocate memory to store information on cells.
-*/
-int buildLandscape(t_Landscape * pLandscape, t_Params * pParams)
-{
- int bRet;
-
- /* blank everything out */
- fprintf(stdout, "entered buildLandscape()\n");
- bRet = 0;
-
- pLandscape->asCells = NULL;
- pLandscape->asUndev = NULL;
- pLandscape->undevSites = 0;
- pLandscape->aParcelSizes = NULL;
- pLandscape->parcelSizes = 0;
-
- /* set initial values across the landscape */
-
- pLandscape->maxX = pParams->xSize;
- pLandscape->maxY = pParams->ySize;
- pLandscape->totalCells = pLandscape->maxX * pLandscape->maxY;
- pLandscape->asCells =
- (t_Cell *) malloc(sizeof(t_Cell) * pLandscape->totalCells);
- if (pLandscape->asCells) {
- if (pParams->num_Regions > 1) {
- /* just allocate enough space for every cell to be undev */
- pLandscape->asUndev = (t_Undev *) malloc(sizeof(t_Undev) * pLandscape->totalCells);
- if (pLandscape->asUndev) {
- bRet = 1;
- }
- }
- else
- bRet = 1;
- }
- return bRet;
-}
-
-void readData4AdditionalVariables(t_Landscape * pLandscape,
- t_Params * pParams)
-{
- int i;
- int ii;
- double val;
-
- for (i = 0; i < pParams->numAddVariables; i++) {
- cout << "reading additional variables File: " <<
- pParams->addVariableFile[i] << "...";
-
- /* open raster map */
- int fd = Rast_open_old(pParams->addVariableFile[i], "");
- RASTER_MAP_TYPE data_type = Rast_get_map_type(fd);
- void *buffer = Rast_allocate_buf(data_type);
-
- ii = 0;
- for (int row = 0; row < pParams->xSize; row++) {
- Rast_get_row(fd, buffer, row, data_type);
- void *ptr = buffer;
-
- for (int col = 0; col < pParams->ySize; col++,
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(data_type))) {
- if (data_type == DCELL_TYPE)
- val = *(DCELL *) ptr;
- else if (data_type == FCELL_TYPE)
- val = *(FCELL *) ptr;
- else {
- val = *(CELL *) ptr;
- }
- if (pLandscape->asCells[ii].nCellType == _CELL_VALID)
- pLandscape->asCells[ii].additionVariable[i] = val;
- else
- pLandscape->asCells[ii].additionVariable[i] = 0.0;
- ii++;
- }
- }
- Rast_close(fd);
- cout << "done" << endl;
- }
-}
-
-void readIndexData(t_Landscape * pLandscape, t_Params * pParams)
-{
- int ii;
- int val;
-
- ii = 0;
- int fd = Rast_open_old(pParams->indexFile, "");
-
- /* TODO: data type should always be CELL */
- RASTER_MAP_TYPE data_type = Rast_get_map_type(fd);
- void *buffer = Rast_allocate_buf(data_type);
-
- cout << "reading index File: " << pParams->indexFile << "...";
- for (int row = 0; row < pParams->xSize; row++) {
- Rast_get_row(fd, buffer, row, data_type);
- void *ptr = buffer;
-
- for (int col = 0; col < pParams->ySize; col++,
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(data_type))) {
- if (Rast_is_null_value(ptr, data_type))
- val = _GIS_NO_DATA_INT; // assign FUTURES representation of NULL value
- else {
- val = *(CELL *) ptr;
- if (val > pParams->num_Regions)
- G_fatal_error(_("Region ID (%d) inconsistent with number of regions (%d) in map <%s>"),
- val, pParams->num_Regions,
- pParams->indexFile);
- if (val < 1)
- G_fatal_error(_("Region ID (%d) is smaller then 1 in map <%s>"),
- val, pParams->indexFile);
- }
-
- pLandscape->asCells[ii].index_region = val;
- ii++;
- }
- }
- cout << "done" << endl;
-}
-
-/**
- Read data from GIS rasters and put in correct places in memory.
-*/
-int readData(t_Landscape * pLandscape, t_Params * pParams)
-{
- // TODO: this function should be rewritten to be nicer (unless storing in memorey will change completely)
- // TODO: fix null handling somewhere
- char *szBuff;
- char szFName[_N_MAX_FILENAME_LEN];
- int bRet, x, y, i, j;
- double dVal;
-
- fprintf(stdout, "entered readData()\n");
- bRet = 0;
- szBuff = (char *)malloc(_N_MAX_DYNAMIC_BUFF_LEN * sizeof(char));
- if (szBuff) {
- for (j = 0; j < 6; j++) {
- switch (j) { /* get correct filename */
- case 0:
- strcpy(szFName, pParams->undevelopedFile);
- break;
- case 1:
- strcpy(szFName, pParams->employAttractionFile);
- break;
- case 2:
- strcpy(szFName, pParams->interchangeDistanceFile);
- break;
- case 3:
- strcpy(szFName, pParams->roadDensityFile);
- break;
- case 4:
- strcpy(szFName, pParams->devPressureFile);
- break;
- case 5:
- strcpy(szFName, pParams->consWeightFile);
- break;
- default:
- fprintf(stderr, "readData(): shouldn't get here...\n");
- break;
- }
- fprintf(stdout, "\t%s...", szFName);
- /* open raster map */
- int fd = Rast_open_old(szFName, "");
-
- i = 0;
-
- x = 0;
- bRet = 1; /* can only get worse from here on in */
- RASTER_MAP_TYPE data_type = Rast_get_map_type(fd);
- void *buffer = Rast_allocate_buf(data_type);
-
- for (int row = 0; row < pParams->xSize; row++) {
- y = 0;
- Rast_get_row(fd, buffer, row, data_type);
- void *ptr = buffer;
-
- for (int col = 0; col < pParams->ySize; col++,
- ptr = G_incr_void_ptr(ptr, Rast_cell_size(data_type))) {
- CELL iVal;
-
- if (data_type == DCELL_TYPE)
- dVal = *(DCELL *) ptr;
- else if (data_type == FCELL_TYPE)
- dVal = *(FCELL *) ptr;
- else {
- iVal = *(CELL *) ptr; // integer value for some rasters
- dVal = iVal;
- }
- if (j == 0) {
- /* check for NO_DATA */
- if (Rast_is_null_value(ptr, data_type)) {
- pLandscape->asCells[i].nCellType =
- _CELL_OUT_OF_COUNTY;
- }
- else {
- pLandscape->asCells[i].nCellType = _CELL_VALID;
- }
- /* and put in the correct place */
- pLandscape->asCells[i].thisX = x;
- pLandscape->asCells[i].thisY = y;
- /* and set the time of development */
- pLandscape->asCells[i].tDeveloped = _GIS_NO_DATA_INT;
- }
- else {
- /* clean up missing data */
- if (Rast_is_null_value(ptr, data_type)) {
- dVal = 0.0;
- }
- }
- // TODO: how do we know that this is set?
- if (pLandscape->asCells[i].nCellType == _CELL_VALID) {
- /* put data into a correct place */
- switch (j) {
- case 0:
- pLandscape->asCells[i].bUndeveloped = iVal;
- pLandscape->asCells[i].bUntouched = 1;
- if (pLandscape->asCells[i].bUndeveloped == 1) {
- pLandscape->asCells[i].tDeveloped =
- _N_NOT_YET_DEVELOPED;
- }
- else {
- /* already developed at the start */
- pLandscape->asCells[i].tDeveloped = 0;
- }
- break;
- case 1:
- pLandscape->asCells[i].employAttraction = dVal;
- break;
- case 2:
- pLandscape->asCells[i].interchangeDistance = dVal;
- break;
- case 3:
- pLandscape->asCells[i].roadDensity = dVal;
- break;
- case 4:
- pLandscape->asCells[i].devPressure = (int)dVal;
- break;
- case 5:
- pLandscape->asCells[i].consWeight = dVal;
- break;
- default:
- fprintf(stderr,
- "readData(): shouldn't get here...\n");
- break;
- }
- }
- y++;
- i++;
- }
-
- x++;
- }
- if (x == pLandscape->maxX) {
-
- }
- else {
- if (bRet) {
- fprintf(stderr, "readData(): x too small\n");
- bRet = 0;
- }
- }
-
- if (y != pLandscape->maxY) {
- if (bRet) {
- G_message(_("x=%d y=%d"), x, y);
- G_message(_(" xSize=%d ySize=%d"), pParams->xSize,
- pParams->ySize);
- G_message(_(" maxX=%d minX=%d"), pLandscape->maxX,
- pLandscape->maxY);
- fprintf(stderr, "readData(): y too small\n");
- bRet = 0;
- }
- }
- else {
- if (bRet && i == pLandscape->totalCells) {
- fprintf(stdout, "done\n");;
- }
- else {
-
- if (bRet) {
- G_message(_("x=%d y=%d"), x, y);
- G_message(_(" xSize=%d ySize=%d"), pParams->xSize,
- pParams->ySize);
- G_message(_(" maxX=%d minX=%d"), pLandscape->maxX,
- pLandscape->maxY);
- fprintf(stderr,
- "readData(): not read in enough points\n");
- bRet = 0;
- }
-
- }
- }
- Rast_close(fd);
-
- }
- free(szBuff);
- }
-
- return bRet;
-}
-
-
-/**
- Compare two t_Undev structures.
-
- To be used in qsort().
-
- \note This sort function is in ascending order
- (i.e. opposite to what we normally want).
-
- \see undevCmpReverse()
-*/
-static int undevCmp(const void *a, const void *b)
-{
- double da = ((t_Undev *) a)->logitVal;
- double db = ((t_Undev *) b)->logitVal;
-
- /* make sure untouched cells preferentially go to top of sort when there is asymmetry between a and b */
- if (((t_Undev *) a)->bUntouched && !((t_Undev *) b)->bUntouched)
- return -1;
- if (((t_Undev *) b)->bUntouched && !((t_Undev *) a)->bUntouched)
- return 1;
- /* then sort by value */
- if (da > db)
- return 1;
- if (da < db)
- return -1;
- return 0;
-}
-
-/**
- Compare two t_Undev structures.
-
- This is the correct way around sorting with undevCmp() and qsort().
- */
-static int undevCmpReverse(const void *a, const void *b)
-{
- return -1 * undevCmp(a, b);
-}
-
-
-/**
- Write current state of developed areas.
-
- Called at end and dumps tDeveloped for all valid cells (NULL for all others)
-
- \param undevelopedAsNull Represent undeveloped areas as NULLs instead of -1
- \param developmentAsOne Represent all developed areas as 1 instead of number
- representing the step when are was developed
- */
-// TODO: add timestamp to maps
-void outputDevRasterStep(t_Landscape * pLandscape, t_Params * pParams,
- const char *rasterNameBase, bool undevelopedAsNull,
- bool developmentAsOne)
-{
- int out_fd = Rast_open_new(rasterNameBase, CELL_TYPE);
- CELL *out_row = Rast_allocate_c_buf();
-
- for (int x = 0; x < pLandscape->maxX; x++) {
- Rast_set_c_null_value(out_row, pLandscape->maxY);
- for (int y = 0; y < pLandscape->maxY; y++) {
- // we don't check return values since we use landscape directly
- int cellId = posFromXY(x, y, pLandscape);
-
- // this handles NULLs on input
- if (pLandscape->asCells[cellId].nCellType == _CELL_VALID) {
- CELL value = pLandscape->asCells[cellId].tDeveloped;
-
- // this handles undeveloped cells
- if (undevelopedAsNull && value == -1)
- continue;
- // this handles developed cells
- if (developmentAsOne)
- value = 1;
- out_row[y] = value;
- }
- }
- Rast_put_c_row(out_fd, out_row);
- }
- Rast_close(out_fd);
-
- struct Colors colors;
-
- Rast_init_colors(&colors);
- CELL val1 = 0;
- // TODO: the map max is 36 for 36 steps, it is correct?
- CELL val2 = pParams->nSteps;
-
- if (developmentAsOne) {
- val1 = 1;
- val2 = 1;
- Rast_add_c_color_rule(&val1, 255, 100, 50, &val2, 255, 100, 50,
- &colors);
- }
- else {
- Rast_add_c_color_rule(&val1, 255, 100, 50, &val2, 255, 255, 0,
- &colors);
- }
- if (!undevelopedAsNull) {
- val1 = -1;
- val2 = -1;
- Rast_add_c_color_rule(&val1, 100, 255, 100, &val2, 100, 255, 100,
- &colors);
- }
-
- const char *mapset = G_find_file2("cell", rasterNameBase, "");
-
- if (mapset == NULL)
- G_fatal_error(_("Raster map <%s> not found"), rasterNameBase);
-
- Rast_write_colors(rasterNameBase, mapset, &colors);
- Rast_free_colors(&colors);
-
- struct History hist;
-
- Rast_short_history(rasterNameBase, "raster", &hist);
- Rast_command_history(&hist);
- // TODO: store also random seed value (need to get it here, global? in Params?)
- Rast_write_history(rasterNameBase, &hist);
-
- G_message(_("Raster map <%s> created"), rasterNameBase);
-}
-
-char *mapNameForStep(const char *basename, const int step, const int maxSteps)
-{
- int digits = log10(maxSteps) + 1;
-
- return G_generate_basename(basename, step, digits, 0);
-}
-
-/**
- Work out how many cells to transition based on the control file.
-*/
-int parseControlLine(char *szBuff, char *szStepLabel, int *pnToConvert)
-{
- char *pPtr;
- int i;
- int bRet;
-
- bRet = 0;
- /* strip newline */
- pPtr = strpbrk(szBuff, "\r\n");
- if (pPtr) {
- pPtr[0] = '\0';
- }
- *pnToConvert = 0;
- i = 0;
- pPtr = strtok(szBuff, " \t");
- while (pPtr && i < 2) {
- switch (i) {
- case 0:
- strcpy(szStepLabel, pPtr);
- break;
- case 1:
- *pnToConvert = atoi(pPtr);
- break;
- }
- pPtr = strtok(NULL, " \t");
- i++;
- }
- /* only process this line of the control file if it is well formed */
- if (i == 2) {
- bRet = 1;
- }
- return bRet;
-}
-
-/**
- Calculate probability of development.
-*/
-
-double getDevProbability(t_Cell * pThis, t_Params * pParams)
-{
- double probAdd;
- int i;
- int id = pThis->index_region;
-
- if (id == -9999)
- return 0;
- id = id - 1;
- probAdd = pParams->dIntercepts[id];
- //cout<<"intercept\t"<<probAdd<<endl;
- probAdd += pParams->dV1[id] * pThis->employAttraction;
- //cout<<"employAttraction: "<<pParams->dV1[id]<<endl;
- //cout<<"employAttraction: "<<pThis->employAttraction<<endl;
- //cout<<"employAttraction: "<<pParams->dV1[id] * pThis->employAttraction<<endl;
- probAdd += pParams->dV2[id] * pThis->interchangeDistance;
- //cout<<"interchangeDistance: "<<pParams->dV2[id]<<endl;
- //cout<<"interchangeDistance: "<<pThis->interchangeDistance<<endl;
- //cout<<"interchangeDistance: "<<pParams->dV2[id] * pThis->interchangeDistance<<endl;
- probAdd += pParams->dV3[id] * pThis->roadDensity;
- //cout<<"roadDensity: "<<pParams->dV3[id]<<endl;
- //cout<<"roadDensity: "<<pThis->roadDensity<<endl;
- //cout<<"roadDensity: "<<pParams->dV3[id] * pThis->roadDensity<<endl;
- probAdd += pParams->dV4[id] * pThis->devPressure;
- //cout<<"devPressure: "<<pParams->dV4[id]<<endl;
- //cout<<"devPressure: "<<pThis->devPressure<<endl;
- //cout<<"devPressure: "<<pParams->dV4[id] * pThis->devPressure<<endl;
- for (i = 0; i < pParams->numAddVariables; i++) {
- probAdd += pParams->addParameters[i][id] * pThis->additionVariable[i];
- //cout<<"additionVariable: "<<i<<"\t"<<pParams->addParameters[i][id]<<endl;
- //cout<<"additionVariable: "<<i<<"\t"<<pThis->additionVariable[i]<<endl;
- //cout<<"additionVariable: "<<i<<"\t"<<pParams->addParameters[i][id] * pThis->additionVariable[i]<<endl;
- }
- probAdd = 1.0 / (1.0 + exp(-probAdd));
- //cout<<"The probability:";
- //cout<<probAdd<<endl;
- return probAdd;
-}
-
-/*
- double getDevProbability(t_Cell *pThis, t_Params *pParams){
- float probAdd; int i;
- int id=pThis->index_region;
- if(id==-9999)
- return 0;
- id=id-1;
- probAdd = pParams->dIntercepts[id];
- probAdd += pParams->dV1[id] * pThis->employAttraction;
- probAdd += pParams->dV2[id] * pThis->interchangeDistance;
- probAdd += pParams->dV3[id] * pThis->roadDensity;
- probAdd += pParams->dV4[id] * pThis->devPressure;
- for(i=0;i<2;i++){
- probAdd+=pParams->addParameters[i][id]*pThis->additionVariable[i];
- }
- if (0<pThis->additionVariable[2]<=400)
- {
- probAdd = probAdd * 0.0005 ;
- }
- else {
- if(400<pThis->additionVariable[2]<=600)
- {probAdd = probAdd * 0.024;}
- else {
- if(600<pThis->additionVariable[2]<=800)
- {probAdd = probAdd * 0.912 ;}
- else{
- if(800<pThis->additionVariable[2]<=1000)
- {probAdd = probAdd * 0.055 ;}
- else{
- if(1000<pThis->additionVariable[2]<=1200)
- {probAdd = probAdd * 0.006 ;}
- else{
- if(1200<pThis->additionVariable[2]<=1400)
- {probAdd = probAdd * 0.002 ;}
- else{
-
- {probAdd = probAdd * 0.0005;}
-
- }
- }
- }
-
- }
-
- }
-
-
- }
-
-
- probAdd = 1.0/(1.0 + exp(-probAdd));
- return probAdd;
- }
-
- */
-
-
-/**
- Recalculate probabilities for each unconverted cell.
-
- Called each timestep.
-*/
-void findAndSortProbs(t_Landscape * pLandscape, t_Params * pParams,
- int nToConvert)
-{
- int i, lookupPos;
- t_Cell *pThis;
-
- /* update calcs */
- fprintf(stdout, "\t\trecalculating probabilities\n");
- pLandscape->undevSites = 0;
- for (i = 0; i < pLandscape->totalCells; i++) {
- pThis = &(pLandscape->asCells[i]);
- if (pThis->nCellType == _CELL_VALID) {
- if (pThis->bUndeveloped) {
- if (pThis->consWeight > 0.0) {
- /* note that are no longer just storing the logit value, but instead the probability (allows consWeight to affect sort order) */
- pLandscape->asUndev[pLandscape->undevSites].cellID = i;
- pLandscape->asUndev[pLandscape->undevSites].logitVal =
- getDevProbability(pThis, pParams);
- /* lookup table of probabilities is applied before consWeight */
- if (pParams->nAlgorithm == _N_ALGORITHM_STOCHASTIC_II) {
- /* replace with value from lookup table */
- lookupPos =
- (int)(pLandscape->asUndev[pLandscape->undevSites].
- logitVal * (pParams->nProbLookup - 1));
- pLandscape->asUndev[pLandscape->undevSites].logitVal =
- pParams->adProbLookup[lookupPos];
- // fprintf(stdout, "%f %d %f\n", pLandscape->asUndev[pLandscape->undevSites].logitVal, lookupPos, pParams->adProbLookup[lookupPos]);
- }
- // multiply consweight
- pLandscape->asUndev[pLandscape->undevSites].logitVal *=
- pThis->consWeight;
- pLandscape->asUndev[pLandscape->undevSites].bUntouched = pThis->bUntouched; /* need to store this to put correct elements near top of list */
- if (pLandscape->asUndev[pLandscape->undevSites].logitVal >
- 0.0) {
- /* only add one more to the list if have a positive probability */
- pLandscape->undevSites++;
- }
- }
- }
- }
- }
- /* downweight the devPressure if necessary (do not do in first step) */
- /* doing it here means that last time step values have full weight */
- if (pParams->nAlgorithm == _N_ALGORITHM_STOCHASTIC_I) {
- if (pParams->dDevPersistence < 1.0) {
- fprintf(stdout, "\t\tdownweighting development pressure\n");
-
- for (i = 0; i < pLandscape->totalCells; i++) {
- pThis = &(pLandscape->asCells[i]);
- if (pThis->nCellType == _CELL_VALID) {
- /* only need to bother downweighting on cells that can still convert */
- if (pThis->bUndeveloped) {
- pThis->devPressure =
- (int)((double)pThis->devPressure *
- pParams->dDevPersistence);
- }
- }
- }
- }
- }
- /* sort */
- /* can only be zero for algorithm=stochastic_ii */
- if (pParams->sortProbs) {
- fprintf(stdout, "\t\tsorting %d unconserved undeveloped sites\n",
- pLandscape->undevSites);
- qsort(pLandscape->asUndev, pLandscape->undevSites, sizeof(t_Undev),
- undevCmpReverse);
- }
- else {
- fprintf(stdout, "\t\tskipping sort as choosing cells randomly\n");
- }
- //calculate cumulative probability
- // From Wenwu Tang
- double sum = pLandscape->asUndev[0].logitVal;
-
- for (i = 1; i < pLandscape->undevSites; i++) {
- pLandscape->asUndev[i].cumulProb =
- pLandscape->asUndev[i - 1].cumulProb +
- pLandscape->asUndev[i].logitVal;
- sum = sum + pLandscape->asUndev[i].logitVal;
- }
- for (i = 0; i < pLandscape->undevSites; i++) {
- pLandscape->asUndev[i].cumulProb =
- pLandscape->asUndev[i].cumulProb / sum;
- }
-}
-
-/**
- Create patches based on spiraling around.
-
- \deprecated
-*/
-int fillValidNeighbourList(int nThisID, t_Landscape * pLandscape,
- t_Params * pParams, int *anToConvert,
- int nWantToConvert, int bAllowTouched,
- int bDeterministic)
-{
- int nTried, nFound, x, y, upDown, stopAt, countMove, stepMove, nPos,
- nSkipped;
-
- anToConvert[0] = nThisID;
- nSkipped = 0;
- nFound = 1;
- nTried = 1;
- x = pLandscape->asCells[nThisID].thisX;
- y = pLandscape->asCells[nThisID].thisY;
- stopAt = 0;
- upDown = 0;
- stepMove = -1;
- while (nFound < nWantToConvert &&
- ((double)nWantToConvert * pParams->giveUpRatio) > nTried) {
- countMove = 0;
- upDown = !upDown;
- if (upDown) {
- stopAt++;
- stepMove *= -1;
- }
- while (countMove < stopAt && nFound < nWantToConvert &&
- ((double)nWantToConvert * pParams->giveUpRatio) > nTried) {
- if (upDown) {
- x += stepMove;
- }
- else {
- y += stepMove;
- }
- // fprintf(stdout, "%d %d\n", x, y);
- nPos = posFromXY(x, y, pLandscape);
- if (nPos != _CELL_OUT_OF_RANGE) {
- if (pLandscape->asCells[nPos].nCellType == _CELL_VALID) {
- if (pLandscape->asCells[nPos].bUndeveloped) {
- if (bAllowTouched ||
- pLandscape->asCells[nPos].bUntouched) {
- if (bDeterministic ||
- (uniformRandom() < pParams->dProbWeight)) {
- if (pLandscape->asCells[nPos].consWeight >
- 0.0) {
- anToConvert[nFound] = nPos;
- nFound++;
- }
- }
- else {
- pLandscape->asCells[nPos].bUntouched = 0;
- nSkipped++;
- }
- }
- }
- }
- }
- nTried++;
- countMove++;
- }
- }
- // fprintf(stdout, "\tskipped=%d\n", nSkipped);
- return nFound;
-}
-
-/**
- Helper structur for neighbour grabbing algorithm
- */
-typedef struct
-{
- double probAdd;
- int cellID;
- int newInList;
- // Wenwu Tang
- /// distance to the center
- double distance;
-} t_candidateNeighbour;
-
-/**
- Helper structur for neighbour grabbing algorithm
-*/
-typedef struct
-{
- double maxProb;
- int nCandidates;
- int nSpace;
- t_candidateNeighbour *aCandidates;
-} t_neighbourList;
-
-#define _N_NEIGHBOUR_LIST_BLOCK_SIZE 20
-
-int addNeighbourIfPoss(int x, int y, t_Landscape * pLandscape,
- t_neighbourList * pNeighbours, t_Params * pParams)
-{
- int i, thisPos, mustAdd, listChanged, lookupPos;
- double probAdd;
- t_Cell *pThis;
-
- listChanged = 0;
- thisPos = posFromXY(x, y, pLandscape);
- if (thisPos != _CELL_OUT_OF_RANGE) {
- pThis = &(pLandscape->asCells[thisPos]);
- if (pThis->nCellType == _CELL_VALID) {
- pThis->bUntouched = 0;
- if (pThis->bUndeveloped) {
- if (pThis->consWeight > 0.0) {
- /* need to add this cell... */
-
- /* ...either refresh its element in list if already there */
- mustAdd = 1;
- for (i = 0; mustAdd && i < pNeighbours->nCandidates; i++) {
- if (pNeighbours->aCandidates[i].cellID == thisPos) {
- pNeighbours->aCandidates[i].newInList = 1;
- mustAdd = 0;
- listChanged = 1;
- }
- }
- /* or add it on the end, allocating space if necessary */
- if (mustAdd) {
- if (pNeighbours->nCandidates == pNeighbours->nSpace) {
- pNeighbours->nSpace +=
- _N_NEIGHBOUR_LIST_BLOCK_SIZE;
- pNeighbours->aCandidates =
- (t_candidateNeighbour *)
- realloc(pNeighbours->aCandidates,
- pNeighbours->nSpace *
- sizeof(t_candidateNeighbour));
- if (!pNeighbours->aCandidates) {
- fprintf(stderr,
- "memory error in addNeighbourIfPoss()\n...exiting\n");
- exit(EXIT_FAILURE);
- }
- }
- pNeighbours->aCandidates[pNeighbours->nCandidates].
- cellID = thisPos;
- pNeighbours->aCandidates[pNeighbours->nCandidates].
- newInList = 1;
-
- /* note duplication of effort in here recalculating the probabilities, but didn't store them in an accessible way */
- /*
- probAdd = pParams->dIntercept;
- probAdd += pParams->dDevPressure * pThis->devPressure;
- probAdd += pParams->dEmployAttraction * pThis->employAttraction;
- probAdd += pParams->dInterchangeDistance * pThis->interchangeDistance;
- probAdd += pParams->dRoadDensity * pThis->roadDensity;
- probAdd = 1.0/(1.0 + exp(probAdd));
- */
- probAdd = getDevProbability(pThis, pParams);
- /* replace with value from lookup table */
- lookupPos =
- (int)(probAdd * (pParams->nProbLookup - 1));
- probAdd = pParams->adProbLookup[lookupPos];
- probAdd *= pThis->consWeight;
- /* multiply by the patch factor to help control shapes of patches */
- // pF now is set to 1, so it won't affect the result.or this can be deleted // WT
- probAdd *= pParams->patchFactor;
- pNeighbours->aCandidates[pNeighbours->nCandidates].
- probAdd = probAdd;
- /* only actually add it if will ever transition */
- if (probAdd > 0.0) {
- pNeighbours->nCandidates++;
- if (probAdd > pNeighbours->maxProb) {
- pNeighbours->maxProb = probAdd;
- }
- listChanged = 1;
- }
- }
- }
- }
- }
- }
- return listChanged;
-}
-
-// From Wenwu
-
-/**
- Sort according to the new flag and conversion probability.
-*/
-static int sortNeighbours(const void *pOne, const void *pTwo)
-{
- t_candidateNeighbour *pNOne = (t_candidateNeighbour *) pOne;
- t_candidateNeighbour *pNTwo = (t_candidateNeighbour *) pTwo;
-
- /*
- if(pNOne->newInList > pNTwo->newInList)
- {
- return -1;
- }
- if(pNTwo->newInList > pNOne->newInList)
- {
- return 1;
- }
- */
- float p = rand() / (float)RAND_MAX;
-
- /*
- if(p<0.0){
- if(pNOne->distance> pNTwo->distance)
- {
- return 1;
- }
- if(pNTwo->distance > pNOne->distance)
- {
- return -1;
- }
- }
- */
- float alpha = 0.5; // assign this from parameter file
-
- alpha = (sParams.patchMean) - (sParams.patchRange) * 0.5;
- alpha += p * sParams.patchRange;
- //alpha=0;
- float p1, p2;
-
- p1 = pNOne->probAdd / pow(pNOne->distance, alpha);
- p2 = pNTwo->probAdd / pow(pNTwo->distance, alpha);
-
- //p1=alpha*pNOne->probAdd+(1-alpha)*pow(pNOne->distance,-2);
- //p2=alpha*pNTwo->probAdd+(1-alpha)*pow(pNTwo->distance,-2);
-
- if (p1 > p2) {
- return -1;
- }
- if (p2 > p1) {
- return 1;
- }
- return 0;
-}
-
-/**
- randomly shuffle the 4 neighbors
-*/
-void Shuffle4(int *mat)
-{
- int proba[4], flag[4];
- int size = 4, i;
-
- for (i = 0; i < size; i++) {
- proba[i] = rand() % 1000;
- flag[i] = 1;
- }
-
- int numChecked = 0;
- int max;
- int index = 0;
-
- while (numChecked != 4) {
- max = -9999;
- for (i = 0; i < size; i++) {
- if (flag[i] != 0) {
- if (proba[i] > max) {
- max = proba[i];
- index = i;
- }
- }
- }
- flag[index] = 0;
- mat[numChecked] = index;
- numChecked++;
- }
-}
-
-void findNeighbours(int nThisID, t_Landscape * pLandscape,
- t_neighbourList * pNeighbours, t_Params * pParams)
-{
- t_Cell *pThis;
- int listChanged = 0;
-
- // int idmat[4]; idmat[0]=0;idmat[1]=1;idmat[2]=2;idmat[3]=3;
- /* try to add the neighbours of this one */
- pThis = &(pLandscape->asCells[nThisID]);
- //note: if sorted, then shuffle is no need
- /*
- Shuffle4(&idmat[0]);
- for(i=0;i<4;i++){
- if(idmat[i]==0)
- listChanged+=addNeighbourIfPoss(pThis->thisX-1,pThis->thisY,pLandscape,pNeighbours,pParams);
- else if(idmat[i]==1)
- listChanged+=addNeighbourIfPoss(pThis->thisX+1,pThis->thisY,pLandscape,pNeighbours,pParams);
- else if(idmat[i]==2)
- listChanged+=addNeighbourIfPoss(pThis->thisX,pThis->thisY-1,pLandscape,pNeighbours,pParams);
- else if (idmat[i]==3)
- listChanged+=addNeighbourIfPoss(pThis->thisX,pThis->thisY+1,pLandscape,pNeighbours,pParams);
- }
- */
-
- listChanged += addNeighbourIfPoss(pThis->thisX - 1, pThis->thisY, pLandscape, pNeighbours, pParams); // left
- listChanged += addNeighbourIfPoss(pThis->thisX + 1, pThis->thisY, pLandscape, pNeighbours, pParams); // right
- listChanged += addNeighbourIfPoss(pThis->thisX, pThis->thisY - 1, pLandscape, pNeighbours, pParams); // down
- listChanged += addNeighbourIfPoss(pThis->thisX, pThis->thisY + 1, pLandscape, pNeighbours, pParams); // up
- if (sParams.numNeighbors == 8) {
- listChanged +=
- addNeighbourIfPoss(pThis->thisX - 1, pThis->thisY - 1, pLandscape,
- pNeighbours, pParams);
- listChanged +=
- addNeighbourIfPoss(pThis->thisX - 1, pThis->thisY + 1, pLandscape,
- pNeighbours, pParams);
- listChanged +=
- addNeighbourIfPoss(pThis->thisX + 1, pThis->thisY - 1, pLandscape,
- pNeighbours, pParams);
- listChanged +=
- addNeighbourIfPoss(pThis->thisX + 1, pThis->thisY + 1, pLandscape,
- pNeighbours, pParams);
- }
- /* if anything has been altered then resort */
-
- if (listChanged) {
- //qsort(pNeighbours->aCandidates,pNeighbours->nCandidates,sizeof(t_candidateNeighbour),sortNeighbours);
- }
-}
-
-double getDistance(int id1, int id2, t_Landscape * pLandscape)
-{
- double x1, x2, y1, y2, result = 0;
-
- x1 = pLandscape->asCells[id1].thisX;
- x2 = pLandscape->asCells[id2].thisX;
- y1 = pLandscape->asCells[id1].thisY;
- y2 = pLandscape->asCells[id2].thisY;
- result = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
- return result;
-}
-
-double getDistance1(double x1, double y1, double x2, double y2)
-{
- double result = 0;
-
- result = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
- return result;
-}
-
-int newPatchFinder(int nThisID, t_Landscape * pLandscape, t_Params * pParams,
- int *anToConvert, int nWantToConvert)
-{
- int bTrying, i, nFound, j;
- t_neighbourList sNeighbours;
-
- memset(&sNeighbours, 0, sizeof(t_neighbourList));
- anToConvert[0] = nThisID;
- pLandscape->asCells[nThisID].bUntouched = 0;
- pLandscape->asCells[nThisID].bUndeveloped = 0;
- nFound = 1;
- findNeighbours(nThisID, pLandscape, &sNeighbours, pParams);
- for (i = 0; i < sNeighbours.nCandidates; i++) {
- sNeighbours.aCandidates[i].distance =
- getDistance(nThisID, sNeighbours.aCandidates[i].cellID,
- pLandscape);
- }
-
- while (nFound < nWantToConvert && sNeighbours.nCandidates > 0) {
- i = 0;
- bTrying = 1;
- while (bTrying) {
- if (uniformRandom() < sNeighbours.aCandidates[i].probAdd) {
- anToConvert[nFound] = sNeighbours.aCandidates[i].cellID;
- /* flag that it is about to develop */
- pLandscape->asCells[anToConvert[nFound]].bUndeveloped = 0;
- /* remove this one from the list by copying down everything above it */
- for (j = i + 1; j < sNeighbours.nCandidates; j++) {
- sNeighbours.aCandidates[j - 1].cellID =
- sNeighbours.aCandidates[j].cellID;
- sNeighbours.aCandidates[j - 1].newInList =
- sNeighbours.aCandidates[j].newInList;
- sNeighbours.aCandidates[j - 1].probAdd =
- sNeighbours.aCandidates[j].probAdd;
- }
- /* reduce the size of the list */
- sNeighbours.nCandidates--;
- sNeighbours.maxProb = 0.0;
- for (j = 0; j < sNeighbours.nCandidates; j++) {
- if (sNeighbours.aCandidates[j].probAdd >
- sNeighbours.maxProb) {
- sNeighbours.maxProb =
- sNeighbours.aCandidates[j].probAdd;
- }
- }
- /* add its neighbours */
- findNeighbours(anToConvert[nFound], pLandscape, &sNeighbours,
- pParams);
- nFound++;
- bTrying = 0;
- }
- else {
- sNeighbours.aCandidates[i].newInList = 0;
- }
- if (bTrying) {
- i++;
- if (i == sNeighbours.nCandidates) {
- i = 0;
- }
- }
- else {
- /* always resort the list if have just added, to let new elements bubble to the top */
- if (sNeighbours.nCandidates > 0) {
- // int z;
- for (i = 0; i < sNeighbours.nCandidates; i++) {
- sNeighbours.aCandidates[i].distance =
- getDistance(nThisID,
- sNeighbours.aCandidates[i].cellID,
- pLandscape);
- }
-
- qsort(sNeighbours.aCandidates, sNeighbours.nCandidates,
- sizeof(t_candidateNeighbour), sortNeighbours);
-#if 0
- for (z = 0; z < sNeighbours.nCandidates; z++) {
- fprintf(stdout, "%d %d %f\n", z,
- sNeighbours.aCandidates[z].newInList,
- sNeighbours.aCandidates[z].probAdd);
- }
- fprintf(stdout, "***\n");
-#endif
- }
- }
- }
- }
- if (sNeighbours.nSpace) { /* free any allocated memory */
- free(sNeighbours.aCandidates);
- }
- // fprintf(stdout, "looking for %d found %d\n", nWantToConvert,nFound);
- return nFound;
-}
-
-int convertCells(t_Landscape * pLandscape, t_Params * pParams, int nThisID,
- int nStep, int bAllowTouched, int bDeterministic)
-{
- int i, x, y, nCell, nDone, nToConvert, nWantToConvert;
- int *anToConvert;
- t_Cell *pThis, *pThat;
-
- nDone = 0;
- /* in here need to build list of near neigbours to loop over */
- nWantToConvert = pLandscape->aParcelSizes[(int)
- (uniformRandom() *
- pLandscape->parcelSizes)];
- anToConvert = (int *)malloc(sizeof(int) * nWantToConvert);
- if (anToConvert) {
- /* in here goes code to fill up list of neighbours */
- if (pParams->nAlgorithm == _N_ALGORITHM_STOCHASTIC_II) {
- nToConvert =
- newPatchFinder(nThisID, pLandscape, pParams, anToConvert,
- nWantToConvert);
- }
- else {
- nToConvert =
- fillValidNeighbourList(nThisID, pLandscape, pParams,
- anToConvert, nWantToConvert,
- bAllowTouched, bDeterministic);
- }
- /* will actually always be the case */
- if (nToConvert > 0) {
- // fprintf(stdout, "wanted %d, found %d\n", nWantToConvert, nToConvert);
- for (i = 0; i < nToConvert; i++) {
- /* convert, updating dev press on neighbours */
- pThis = &(pLandscape->asCells[anToConvert[i]]);
- pThis->bUndeveloped = 0;
- pThis->bUntouched = 0;
- pThis->tDeveloped = nStep + 1;
- nDone++;
- // fprintf(stdout, "%d total=%f dev=%f employment=%f interchange=%f roads=%f\n", pLandscape->asUndev[i].cellID, pLandscape->asUndev[i].logitVal, pParams->dDevPressure * pThis->devPressure, pParams->dEmployAttraction * pThis->employAttraction, pParams->dInterchangeDistance * pThis->interchangeDistance, pParams->dRoadDensity * pThis->roadDensity);
- float dist = 0;
- float force = 0;
-
- for (x = pThis->thisX - pParams->nDevNeighbourhood;
- x <= pThis->thisX + pParams->nDevNeighbourhood; x++) {
- for (y = pThis->thisY - pParams->nDevNeighbourhood;
- y <= pThis->thisY + pParams->nDevNeighbourhood;
- y++) {
- nCell = posFromXY(x, y, pLandscape);
- if (nCell != _CELL_OUT_OF_RANGE) {
- pThat = &(pLandscape->asCells[nCell]);
- dist =
- getDistance1(pThis->thisX, pThis->thisY, x,
- y);
- if (dist > pParams->nDevNeighbourhood)
- continue;
- force = 0;
- if (pParams->devPressureApproach == 1)
- force = 1; // development pressure increases by 1
- else if (pParams->devPressureApproach == 2) {
- if (dist > 0)
- force =
- pParams->scalingFactor / pow(dist,
- pParams->
- alpha);
- }
- else {
- force =
- pParams->scalingFactor * exp(-2 * dist /
- pParams->
- alpha);
- }
- pThat->devPressure = pThat->devPressure + force;
- }
- }
- }
- }
- }
- free(anToConvert);
- }
- return nDone;
-}
-
-/**
- main routine to actually run the model
-*/
-void updateMap(t_Landscape * pLandscape, t_Params * pParams)
-{
- FILE *fIn;
- char *szBuff;
- char szStepLabel[_N_MAX_FILENAME_LEN];
- int i, nToConvert, nStep, nDone, bAllowTouched, nRandTries, nExtra;
- double dProb;
- t_Cell *pThis;
-
- nExtra = 0;
- fprintf(stdout, "entered updateMap()\n");
- fIn = fopen(pParams->controlFile, "rb");
- if (fIn) {
- szBuff = (char *)malloc(_N_MAX_DYNAMIC_BUFF_LEN * sizeof(char));
- if (szBuff) {
- /* start counter at 1 so can distinguish between cells
- which transition on first step and those which already developed */
- nStep = 1;
- while (fgets(szBuff, _N_MAX_DYNAMIC_BUFF_LEN, fIn)) {
- if (parseControlLine(szBuff, szStepLabel, &nToConvert)) {
- fprintf(stdout,
- "\tdoing step %s...controlFile requests conversion of %d cells\n",
- szStepLabel, nToConvert);
- if (nExtra > 0) {
- if (nToConvert - nExtra > 0) {
- nToConvert -= nExtra;
- nExtra = 0;
- }
- else {
- nExtra -= nToConvert;
- nToConvert = 0;
- }
- }
- fprintf(stdout,
- "\t\tafter accounting for extra cells, attempt %d cells\n",
- nToConvert);
- /* if have cells to convert this step */
- if (nToConvert > 0) {
- findAndSortProbs(pLandscape, pParams, nToConvert);
- /* if not enough cells to convert then alter number required */
- if (nToConvert > pLandscape->undevSites) {
- fprintf(stdout,
- "\t\tnot enough undeveloped sites...converting all\n");
- nToConvert = pLandscape->undevSites;
- }
- /* update either in deterministic or stochastic fashion */
- fprintf(stdout, "\t\tupdating map\n");
- switch (pParams->nAlgorithm) {
- case _N_ALGORITHM_DETERMINISTIC: /* deterministic */
- nDone = 0;
- for (i = 0; i < nToConvert && nDone < nToConvert;
- i++) {
- nDone +=
- convertCells(pLandscape, pParams,
- pLandscape->asUndev[i].
- cellID, nStep, 1, 1);
- }
- break;
- case _N_ALGORITHM_STOCHASTIC_I: /* stochastic */
- nDone = 0;
- i = 0;
- bAllowTouched = 0;
- /* loop until done enough cells...might need multiple passes */
- while (nDone < nToConvert) {
- if (i == pLandscape->undevSites) {
- /* if at the end of the grid, just loop around again until done */
- i = 0;
- /* allow previously considered cells if you have to */
- bAllowTouched = 1;
- }
- pThis =
- &(pLandscape->asCells
- [pLandscape->asUndev[i].cellID]);
- if (pThis->bUndeveloped) { /* need to check is still undeveloped */
- if (bAllowTouched || pThis->bUntouched) {
- /* Doug's "back to front" logit */
- // dProb = 1.0/(1.0 + exp(pLandscape->asUndev[i].logitVal));
- dProb =
- pLandscape->asUndev[i].logitVal;
- /* if starting a patch off here */
- if (uniformRandom() <
- pParams->dProbWeight * dProb) {
- nDone +=
- convertCells(pLandscape,
- pParams,
- pLandscape->
- asUndev[i].
- cellID, nStep,
- bAllowTouched,
- 0);
- }
- }
- pThis->bUntouched = 0;
- }
- i++;
- }
- break;
- case _N_ALGORITHM_STOCHASTIC_II: /* stochastic */
- nDone = 0;
- i = 0;
- nRandTries = 0;
- bAllowTouched = 0;
- /* loop until done enough cells...might need multiple passes */
- while (nDone < nToConvert) {
- if (i == pLandscape->undevSites) {
- /* if at the end of the grid, just loop around again until done */
- i = 0;
- /* allow previously considered cells if you have to */
- bAllowTouched = 1;
- }
- /* if tried too many randomly in this step, give up on idea of only letting untouched cells convert */
- if (nRandTries >
- _MAX_RAND_FIND_SEED_FACTOR * nToConvert) {
- bAllowTouched = 1;
- }
- if (pParams->sortProbs) {
- /* if sorted then choose the top cell and do nothing */
- }
- else {
- /* otherwise give a random undeveloped cell a go */
- if (sParams.seedSearch == 1)
- i = (int)(uniformRandom() *
- pLandscape->undevSites);
- // pick one according to their probability
- else
- i = getUnDevIndex(pLandscape);
-
- }
- pThis =
- &(pLandscape->asCells
- [pLandscape->asUndev[i].cellID]);
-
- if (pThis->bUndeveloped) { /* need to check is still undeveloped */
- if (bAllowTouched || pThis->bUntouched) {
- /* Doug's "back to front" logit */
- // dProb = 1.0/(1.0 + exp(pLandscape->asUndev[i].logitVal));
- dProb =
- pLandscape->asUndev[i].logitVal;
- if (uniformRandom() < dProb) {
- nDone +=
- convertCells(pLandscape,
- pParams,
- pLandscape->
- asUndev[i].
- cellID, nStep,
- bAllowTouched,
- 0);
- }
- }
- pThis->bUntouched = 0;
- }
- if (pParams->sortProbs) {
- i++;
- }
- else {
- nRandTries++;
- }
- }
- break;
- default:
- fprintf(stderr, "Unknown algorithm...exiting\n");
- break;
- }
- fprintf(stdout, "\t\tconverted %d sites\n", nDone);
- nExtra += (nDone - nToConvert);
- fprintf(stdout,
- "\t\t%d extra sites knocked off next timestep\n",
- nExtra);
- }
- }
- nStep++; /* next time step */
- }
- /* dump results to a file at the end of the run */
- outputDevRasterStep(pLandscape, pParams, pParams->dumpFile, false,
- false);
- free(szBuff);
- }
- fclose(fIn);
- }
-}
-
-void readDevDemand(t_Params * pParams)
-{
- ifstream f1;
-
- /*
- f1.open(pParams->controlFile);
- int counter=0;
- int val1,val2;
- while(!f1.eof()){
- f1>>val1>>val2;
- pParams->devDemand[counter]=val2;
- counter++;
- }
- pParams->nSteps=counter;
- f1.close();
- */
- //read land demand for each region
- int i;
- int val1;
- char str[200];
-
- f1.open(pParams->controlFileAll);
- f1 >> str >> val1;
- pParams->nSteps = val1;
- f1.getline(str, 200);
- f1.getline(str, 200); //skip the header
- int ii;
-
- for (ii = 0; ii < pParams->nSteps; ii++) {
- f1 >> val1;
- for (i = 0; i < pParams->num_Regions; i++) {
- f1 >> val1;
- pParams->devDemands[i][ii] = val1;
- }
-
- }
- f1.close();
-
-}
-
-void initializeUnDev(t_Landscape * pLandscape, t_Params * pParams)
-{
- int i;
-
- pLandscape->asUndevs.clear();
- pLandscape->asUndevs.resize(pParams->num_Regions,
- std::vector < t_Undev > (MAX_UNDEV_SIZE));
- for (i = 0; i < pParams->num_Regions; i++) {
- pLandscape->num_undevSites[i] = 0;
- }
-}
-
-void finalizeUnDev(t_Landscape * pLandscape, t_Params * pParams)
-{
-
-}
-
-/*
- main routine to run models on multiple regions // WTang
- */
-void updateMapAll(t_Landscape * pLandscape, t_Params * pParams)
-{
- //readDevDemand(pParams);
- int i, j;
-
- //initialize undev arrays
- initializeUnDev(pLandscape, pParams);
-
- //iterate each step (year)
- for (i = 0; i < pParams->nSteps; i++) {
- cout << i << "\t" << pParams->nSteps << endl;
- // for each sub-region, find and update conversion probability (conservation weight applied)
- findAndSortProbsAll(pLandscape, pParams, i);
- for (j = 0; j < pParams->num_Regions; j++)
- //j=1;
- updateMap1(pLandscape, pParams, i, j);
- if (pParams->outputSeries)
- outputDevRasterStep(pLandscape, pParams,
- mapNameForStep(pParams->outputSeries, i,
- pParams->nSteps), true, true);
- }
- outputDevRasterStep(pLandscape, pParams, pParams->dumpFile, false, false);
- finalizeUnDev(pLandscape, pParams);
-}
-
-void updateMap1(t_Landscape * pLandscape, t_Params * pParams, int step,
- int regionID)
-{
- int i, nToConvert, nStep, nDone, bAllowTouched, nRandTries, nExtra;
- double dProb;
- t_Cell *pThis;
-
- nExtra = 0;
-
- nStep = step;
-
- //fprintf(stdout, "entered updateMap()\n");
- nToConvert = pParams->devDemands[regionID][step];
- //fprintf(stdout, "\tdoing step %s...controlFile requests conversion of %d cells\n", szStepLabel, nToConvert);
- if (nExtra > 0) {
- if (nToConvert - nExtra > 0) {
- nToConvert -= nExtra;
- nExtra = 0;
- }
- else {
- nExtra -= nToConvert;
- nToConvert = 0;
- }
- }
- fprintf(stdout,
- "\t\tafter accounting for extra cells, attempt %d cells\n",
- nToConvert);
- /* if have cells to convert this step */
- if (nToConvert > 0) {
- //findAndSortProbs(pLandscape, pParams, nToConvert);
- /* if not enough cells to convert then alter number required */
- if (nToConvert > pLandscape->num_undevSites[regionID]) {
- fprintf(stdout,
- "\t\tnot enough undeveloped sites...converting all\n");
- nToConvert = pLandscape->num_undevSites[regionID];
- }
- /* update either in deterministic or stochastic fashion */
- //fprintf(stdout, "\t\tupdating map\n");
- switch (pParams->nAlgorithm) {
- case _N_ALGORITHM_DETERMINISTIC: /* deterministic */
- nDone = 0;
- for (i = 0; i < nToConvert && nDone < nToConvert; i++) {
- nDone +=
- convertCells(pLandscape, pParams,
- pLandscape->asUndev[i].cellID, nStep, 1, 1);
- }
- break;
- case _N_ALGORITHM_STOCHASTIC_I: /* stochastic */
- nDone = 0;
- i = 0;
- bAllowTouched = 0;
- /* loop until done enough cells...might need multiple passes */
- while (nDone < nToConvert) {
- if (i == pLandscape->undevSites) {
- /* if at the end of the grid, just loop around again until done */
- i = 0;
- /* allow previously considered cells if you have to */
- bAllowTouched = 1;
- }
- pThis = &(pLandscape->asCells[pLandscape->asUndev[i].cellID]);
- /* need to check is still undeveloped */
- if (pThis->bUndeveloped) {
- if (bAllowTouched || pThis->bUntouched) {
- /* Doug's "back to front" logit */
- // dProb = 1.0/(1.0 + exp(pLandscape->asUndev[i].logitVal));
- dProb = pLandscape->asUndev[i].logitVal;
- /* if starting a patch off here */
- if (uniformRandom() < pParams->dProbWeight * dProb) {
- nDone +=
- convertCells(pLandscape, pParams,
- pLandscape->asUndev[i].cellID,
- nStep, bAllowTouched, 0);
- }
- }
- pThis->bUntouched = 0;
- }
- i++;
- }
- break;
- case _N_ALGORITHM_STOCHASTIC_II: /* stochastic */
- nDone = 0;
- i = 0;
- nRandTries = 0;
- bAllowTouched = 0;
- /* loop until done enough cells...might need multiple passes */
- while (nDone < nToConvert) {
- if (i == pLandscape->num_undevSites[regionID]) {
- /* if at the end of the grid, just loop around again until done */
- i = 0;
- /* allow previously considered cells if you have to */
- bAllowTouched = 1;
- }
- /* if tried too many randomly in this step, give up on idea of only letting untouched cells convert */
- if (nRandTries > _MAX_RAND_FIND_SEED_FACTOR * nToConvert) {
- bAllowTouched = 1;
- }
- if (pParams->sortProbs) {
- /* if sorted then choose the top cell and do nothing */
- }
- else {
- /* otherwise give a random undeveloped cell a go */
- if (sParams.seedSearch == 1)
- i = (int)(uniformRandom() *
- pLandscape->num_undevSites[regionID]);
- //pick one according to their probability
- else
- G_verbose_message("nDone=%d, toConvert=%d", nDone,
- nToConvert);
- i = getUnDevIndex1(pLandscape, regionID);
- }
- pThis =
- &(pLandscape->asCells
- [pLandscape->asUndevs[regionID][i].cellID]);
-
- /* need to check is still undeveloped */
- if (pThis->bUndeveloped) {
- if (bAllowTouched || pThis->bUntouched) {
- /* Doug's "back to front" logit */
- // dProb = 1.0/(1.0 + exp(pLandscape->asUndev[i].logitVal));
- dProb = pLandscape->asUndevs[regionID][i].logitVal;
- if (uniformRandom() < dProb) {
- nDone +=
- convertCells(pLandscape, pParams,
- pLandscape->asUndevs[regionID]
- [i].cellID, nStep, bAllowTouched,
- 0);
- }
- }
- pThis->bUntouched = 0;
- }
- if (pParams->sortProbs) {
- i++;
- }
- else {
- nRandTries++;
- }
- }
- break;
- default:
- fprintf(stderr, "Unknown algorithm...exiting\n");
- break;
- }
- fprintf(stdout, "\t\tconverted %d sites\n", nDone);
- nExtra += (nDone - nToConvert);
- fprintf(stdout, "\t\t%d extra sites knocked off next timestep\n",
- nExtra);
- }
-}
-
-
-/**
- Check Doug's calculation of devPressure.
-
- No longer called.
- */
-void testDevPressure(t_Landscape * pLandscape, t_Params * pParams)
-{
- int x, y, xN, yN, cellID, nCell, testPressure;
-
- for (y = 0; y < pLandscape->maxY; y++) {
- for (x = 0; x < pLandscape->maxX; x++) {
- cellID = posFromXY(x, y, pLandscape);
- if (cellID != _CELL_OUT_OF_RANGE) { /* should never happen */
- if (pLandscape->asCells[cellID].nCellType == _CELL_VALID) {
- if (pLandscape->asCells[cellID].bUndeveloped) {
- if (pLandscape->asCells[cellID].devPressure) {
- fprintf(stdout,
- "(x,y)=(%d,%d) cellType=%d bUndev=%d devPress=%f\n",
- x, y,
- pLandscape->asCells[cellID].nCellType,
- pLandscape->asCells[cellID].bUndeveloped,
- pLandscape->asCells[cellID].devPressure);
- testPressure = 0;
- for (xN = x - pParams->nDevNeighbourhood;
- xN <= x + pParams->nDevNeighbourhood; xN++) {
- for (yN = y - pParams->nDevNeighbourhood;
- yN <= y + pParams->nDevNeighbourhood;
- yN++) {
- nCell = posFromXY(xN, yN, pLandscape);
- if (nCell != _CELL_OUT_OF_RANGE) {
- if (pLandscape->asCells[nCell].
- nCellType == _CELL_VALID) {
- if (pLandscape->asCells[nCell].
- bUndeveloped == 0) {
- testPressure++;
- }
- }
- }
- }
- }
- fprintf(stdout, "\t%d\n", testPressure);
- }
- }
- }
- }
- }
- }
-}
-
-int readParcelSizes(t_Landscape * pLandscape, t_Params * pParams)
-{
- FILE *fIn;
- char *szBuff;
- int nMaxParcels;
-
-
- pLandscape->parcelSizes = 0;
-
- fprintf(stdout, "entered readParcelSizes()\n");
- fIn = fopen(pParams->parcelSizeFile, "rb");
- if (fIn) {
- szBuff = (char *)malloc(_N_MAX_DYNAMIC_BUFF_LEN * sizeof(char));
- if (szBuff) {
- /* just scan the file twice */
- nMaxParcels = 0;
- while (fgets(szBuff, _N_MAX_DYNAMIC_BUFF_LEN, fIn)) {
- nMaxParcels++;
- }
- rewind(fIn);
- if (nMaxParcels) {
- pLandscape->aParcelSizes =
- (int *)malloc(sizeof(int) * nMaxParcels);
- if (pLandscape->aParcelSizes) {
- while (fgets(szBuff, _N_MAX_DYNAMIC_BUFF_LEN, fIn)) {
- pLandscape->aParcelSizes[pLandscape->parcelSizes] =
- atoi(szBuff) * pParams->discountFactor;
- if (pLandscape->aParcelSizes[pLandscape->parcelSizes]) {
- pLandscape->parcelSizes++;
- }
- }
- }
- }
- free(szBuff);
- }
- fclose(fIn);
- }
- return pLandscape->parcelSizes;
-}
-
-int main(int argc, char **argv)
-{
-
- struct
- {
- struct Option
- *controlFile, *employAttractionFile, *interchangeDistanceFile,
- *roadDensityFile, *undevelopedFile, *devPressureFile,
- *consWeightFile, *addVariableFiles, *nDevNeighbourhood,
- *devpotParamsFile, *dumpFile, *outputSeries, *algorithm,
- *dProbWeight, *dDevPersistence, *parcelSizeFile, *discountFactor,
- *giveUpRatio, *probLookupFile, *sortProbs, *patchFactor,
- *patchMean, *patchRange, *numNeighbors, *seedSearch,
- *devPressureApproach, *alpha, *scalingFactor, *num_Regions,
- *indexFile, *controlFileAll, *seed;
-
- } opt;
-
- struct
- {
- struct Flag *generateSeed;
- } flg;
-
-
- G_gisinit(argv[0]);
-
- struct GModule *module = G_define_module();
-
- G_add_keyword(_("raster"));
- G_add_keyword(_("patch growing"));
- G_add_keyword(_("urban"));
- G_add_keyword(_("landscape"));
- G_add_keyword(_("modeling"));
- module->label =
- _("Simulates landuse change using FUTure Urban-Regional Environment Simulation (FUTURES).");
- module->description =
- _("Module uses Patch-Growing Algorithm (PGA) to"
- " simulate urban-rural landscape structure development.");
-
- opt.controlFile = G_define_standard_option(G_OPT_F_INPUT);
- opt.controlFile->key = "control_file";
- opt.controlFile->required = NO;
- opt.controlFile->label =
- _("File containing information on how many cells to transition and when");
- opt.controlFile->description =
- _("Needed only for single region/zone run");
-
- opt.undevelopedFile = G_define_standard_option(G_OPT_R_INPUT);
- opt.undevelopedFile->key = "undeveloped";
- opt.undevelopedFile->required = YES;
- opt.undevelopedFile->description =
- _("Files containing the information to read in");
-
- opt.employAttractionFile = G_define_standard_option(G_OPT_R_INPUT);
- opt.employAttractionFile->key = "employ_attraction";
- opt.employAttractionFile->required = YES;
- opt.employAttractionFile->description =
- _("Files containing the information to read in");
-
- opt.interchangeDistanceFile = G_define_standard_option(G_OPT_R_INPUT);
- opt.interchangeDistanceFile->key = "interchange_distance";
- opt.interchangeDistanceFile->required = YES;
- opt.interchangeDistanceFile->description =
- _("Files containing the information to read in");
-
- opt.roadDensityFile = G_define_standard_option(G_OPT_R_INPUT);
- opt.roadDensityFile->key = "road_density";
- opt.roadDensityFile->required = YES;
- opt.roadDensityFile->description =
- _("Files containing the information to read in");
-
- opt.devPressureFile = G_define_standard_option(G_OPT_R_INPUT);
- opt.devPressureFile->key = "development_pressure";
- opt.devPressureFile->required = YES;
- opt.devPressureFile->description =
- _("Files containing the information to read in");
-
- opt.consWeightFile = G_define_standard_option(G_OPT_R_INPUT);
- opt.consWeightFile->key = "cons_weight";
- opt.consWeightFile->required = YES;
- opt.consWeightFile->description =
- _("Files containing the information to read in");
-
- opt.addVariableFiles = G_define_standard_option(G_OPT_R_INPUT);
- opt.addVariableFiles->key = "additional_variable_files";
- opt.addVariableFiles->required = YES;
- opt.addVariableFiles->multiple = YES;
- opt.addVariableFiles->description = _("additional variables");
-
- opt.nDevNeighbourhood = G_define_option();
- opt.nDevNeighbourhood->key = "n_dev_neighbourhood";
- opt.nDevNeighbourhood->type = TYPE_INTEGER;
- opt.nDevNeighbourhood->required = YES;
- opt.nDevNeighbourhood->description =
- _("Size of square used to recalculate development pressure");
-
- opt.devpotParamsFile = G_define_standard_option(G_OPT_F_INPUT);
- opt.devpotParamsFile->key = "devpot_params";
- opt.devpotParamsFile->required = YES;
- opt.devpotParamsFile->label =
- _("Development potential parameters for each region");
- opt.devpotParamsFile->description =
- _("Each line should contain region ID followed"
- " by parameters. Values are separated by whitespace (spaces or tabs)."
- " First line is ignored, so it can be used for header");
-
- opt.algorithm = G_define_option();
- opt.algorithm->key = "algorithm";
- opt.algorithm->type = TYPE_STRING;
- opt.algorithm->required = YES;
- opt.algorithm->options = "deterministic,stochastic1,stochastic2";
- opt.algorithm->description =
- _("Parameters controlling the algorithm to use");
-
- opt.dProbWeight = G_define_option();
- opt.dProbWeight->key = "prob_weight";
- opt.dProbWeight->type = TYPE_DOUBLE;
- opt.dProbWeight->required = NO;
- opt.dProbWeight->label = _("Parameters controlling the algorithm to use");
- opt.dProbWeight->description = _("(only relevant if nAlgorithm=2)"
- " the probabilities are multiplied"
- " by this factor before comparing with random numbers...increases randomness"
- " of the simulation as with probs like 0.95 is more or less deterministic");
- opt.dProbWeight->guisection = _("Stochastic 1");
-
- opt.dDevPersistence = G_define_option();
- opt.dDevPersistence->key = "dev_neighbourhood";
- opt.dDevPersistence->type = TYPE_DOUBLE;
- opt.dDevPersistence->required = NO;
- opt.dDevPersistence->label =
- _("Parameters controlling the algorithm to use");
- opt.dDevPersistence->description =
- _("(only relevant if nAlgorithm=2) the devPressures"
- " are multiplied by this factor on each timestep, meaning older development is"
- " downweighted and more recent development most important...should lead"
- " to clustering of new development");
- opt.dDevPersistence->guisection = _("Stochastic 1");
-
- opt.parcelSizeFile = G_define_standard_option(G_OPT_F_INPUT);
- opt.parcelSizeFile->key = "parcel_size_file";
- opt.parcelSizeFile->required = YES;
- opt.parcelSizeFile->description =
- _("File containing information on the parcel size to use");
-
- opt.discountFactor = G_define_option();
- opt.discountFactor->key = "discount_factor";
- opt.discountFactor->type = TYPE_DOUBLE;
- opt.discountFactor->required = YES;
- opt.discountFactor->description = _("discount factor of patch size");
-
- opt.giveUpRatio = G_define_option();
- opt.giveUpRatio->key = "give_up_ratio";
- opt.giveUpRatio->type = TYPE_DOUBLE;
- opt.giveUpRatio->required = NO;
- opt.giveUpRatio->label = _("Give up ratio");
- opt.giveUpRatio->description = _("(only relevant if nAlgorithm=2) give up"
- " spiralling around when examined this factor of sites too many");
- opt.giveUpRatio->guisection = _("Stochastic 1");
-
- /* stochastic 2 algorithm */
-
- opt.probLookupFile = G_define_standard_option(G_OPT_F_INPUT);
- opt.probLookupFile->key = "probability_lookup_file";
- opt.probLookupFile->required = NO;
- opt.probLookupFile->label =
- _("File containing lookup table for probabilities");
- opt.probLookupFile->description =
- _("Format is tightly constrained. See documentation.");
- opt.probLookupFile->guisection = _("Stochastic 2");
-
- opt.sortProbs = G_define_option();
- opt.sortProbs->key = "sort_probs";
- opt.sortProbs->type = TYPE_INTEGER;
- opt.sortProbs->required = NO;
- opt.sortProbs->description =
- _("Whether or not to sort the list of undeveloped cells before choosing patch seeds");
- opt.sortProbs->guisection = _("Stochastic 2");
-
- opt.patchFactor = G_define_option();
- opt.patchFactor->key = "patch_factor";
- opt.patchFactor->type = TYPE_DOUBLE;
- opt.patchFactor->required = NO;
- opt.patchFactor->description =
- _("when building patches, multiply all probabilities by this"
- " factor (will controls shape of patches to some extent, with higher"
- " numbers more circular and lower numbers more linear)");
- opt.patchFactor->guisection = _("Stochastic 2");
-
- opt.patchMean = G_define_option();
- opt.patchMean->key = "patch_mean";
- opt.patchMean->type = TYPE_DOUBLE;
- opt.patchMean->required = NO;
- opt.patchMean->description =
- _("patch_mean and patch_range are now used to control patch shape");
- opt.patchMean->guisection = _("Stochastic 2");
-
- opt.patchRange = G_define_option();
- opt.patchRange->key = "patch_range";
- opt.patchRange->type = TYPE_DOUBLE;
- opt.patchRange->required = NO;
- opt.patchRange->description =
- _("patch_mean and patch_range are now used to control patch shape");
- opt.patchRange->guisection = _("Stochastic 2");
-
- opt.numNeighbors = G_define_option();
- opt.numNeighbors->key = "num_neighbors";
- opt.numNeighbors->type = TYPE_INTEGER;
- opt.numNeighbors->required = NO;
- opt.numNeighbors->options = "4,8";
- opt.numNeighbors->description =
- _("The number of neighbors to be used for patch generation (4 or 8)");
- opt.numNeighbors->guisection = _("Stochastic 2");
-
- opt.seedSearch = G_define_option();
- opt.seedSearch->key = "seed_search";
- opt.seedSearch->type = TYPE_INTEGER;
- opt.seedSearch->required = NO;
- opt.seedSearch->options = "1,2";
- opt.seedSearch->description =
- _("The way that the location of a seed is determined");
- opt.seedSearch->guisection = _("Stochastic 2");
-
- opt.devPressureApproach = G_define_option();
- opt.devPressureApproach->key = "development_pressure_approach";
- opt.devPressureApproach->type = TYPE_INTEGER;
- opt.devPressureApproach->required = NO;
- opt.devPressureApproach->options = "1,2,3";
- opt.devPressureApproach->description =
- _("approaches to derive development pressure");
- opt.devPressureApproach->guisection = _("Stochastic 2");
-
- opt.alpha = G_define_option();
- opt.alpha->key = "alpha";
- opt.alpha->type = TYPE_DOUBLE;
- opt.alpha->required = NO;
- opt.alpha->description =
- _("Required for development_pressure_approach 1 and 2");
- opt.alpha->guisection = _("Stochastic 2");
-
- opt.scalingFactor = G_define_option();
- opt.scalingFactor->key = "scaling_factor";
- opt.scalingFactor->type = TYPE_DOUBLE;
- opt.scalingFactor->required = NO;
- opt.scalingFactor->description =
- _("Required for development_pressure_approach 2 and 3");
- opt.scalingFactor->guisection = _("Stochastic 2");
-
- opt.num_Regions = G_define_option();
- opt.num_Regions->key = "num_regions";
- opt.num_Regions->type = TYPE_INTEGER;
- opt.num_Regions->required = NO;
- opt.num_Regions->description =
- _("Number of sub-regions (e.g., counties) to be simulated");
- opt.num_Regions->guisection = _("Stochastic 2");
-
- opt.indexFile = G_define_standard_option(G_OPT_R_INPUT);
- opt.indexFile->key = "index_file";
- opt.indexFile->required = NO;
- opt.indexFile->description = _("File for index of sub-regions");
- opt.indexFile->guisection = _("Stochastic 2");
-
- opt.controlFileAll = G_define_standard_option(G_OPT_F_INPUT);
- opt.controlFileAll->key = "control_file_all";
- opt.controlFileAll->required = NO;
- opt.controlFileAll->description =
- _("Control file with number of cells to convert");
- opt.controlFileAll->guisection = _("Stochastic 2");
-
- opt.dumpFile = G_define_standard_option(G_OPT_R_OUTPUT);
- opt.dumpFile->key = "output";
- opt.dumpFile->required = YES;
- opt.dumpFile->description =
- _("State of the development at the end of simulation");
- opt.dumpFile->guisection = _("Output");
-
- opt.outputSeries = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT);
- opt.outputSeries->key = "output_series";
- opt.outputSeries->required = NO;
- opt.outputSeries->label =
- _("State of the development at after each step");
- opt.outputSeries->guisection = _("Output");
- // TODO: add mutually exclusive?
- // TODO: add flags or options to control values in series and final rasters
-
- opt.seed = G_define_option();
- opt.seed->key = "random_seed";
- opt.seed->type = TYPE_INTEGER;
- opt.seed->required = NO;
- opt.seed->label = _("Seed for random number generator");
- opt.seed->description =
- _("The same seed can be used to obtain same results"
- " or random seed can be generated by other means.");
- opt.seed->guisection = _("Random numbers");
-
- flg.generateSeed = G_define_flag();
- flg.generateSeed->key = 's';
- flg.generateSeed->label =
- _("Generate random seed (result is non-deterministic)");
- flg.generateSeed->description =
- _("Automatically generates random seed for random number"
- " generator (use when you don't want to provide the seed option)");
- flg.generateSeed->guisection = _("Random numbers");
-
- // provided XOR generated
- G_option_exclusive(opt.seed, flg.generateSeed, NULL);
- G_option_required(opt.seed, flg.generateSeed, NULL);
-
- if (G_parser(argc, argv))
- exit(EXIT_FAILURE);
-
- long seed_value;
-
- if (flg.generateSeed->answer) {
- seed_value = G_srand48_auto();
- G_message("Generated random seed (-s): %ld", seed_value);
- }
- if (opt.seed->answer) {
- seed_value = atol(opt.seed->answer);
- // this does nothing since we are not using GRASS random function
- G_srand48(seed_value);
- G_message("Read random seed from %s option: %ld",
- opt.seed->key, seed_value);
- }
- // although GRASS random function is initialized
- // the general one must be done separately
- // TODO: replace all by GRASS random number generator?
- srand(seed_value);
-
- // TODO: move this back to local variables
- //t_Params sParams;
- //t_Landscape sLandscape;
-
- /* blank everything out */
- memset(&sParams, 0, sizeof(t_Params));
-
- /* parameters dependednt on region */
- sParams.xSize = Rast_window_rows();
- sParams.ySize = Rast_window_cols();
- G_message(_("Running on %d rows and %d columns (%d cells)"),
- sParams.xSize, sParams.ySize, sParams.xSize * sParams.ySize);
-
- /* set up parameters */
- sParams.controlFile = opt.controlFile->answer;
- sParams.undevelopedFile = opt.undevelopedFile->answer;
- sParams.employAttractionFile = opt.employAttractionFile->answer;
- sParams.interchangeDistanceFile = opt.interchangeDistanceFile->answer;
- sParams.roadDensityFile = opt.roadDensityFile->answer;
- sParams.devPressureFile = opt.devPressureFile->answer;
- sParams.consWeightFile = opt.consWeightFile->answer;
- sParams.numAddVariables = 0;
- char **answer = opt.addVariableFiles->answers;
- size_t num_answers = 0;
-
- while (opt.addVariableFiles->answers[num_answers]) {
- sParams.addVariableFile[num_answers] = *answer;
- num_answers++;
- }
- sParams.numAddVariables = num_answers;
- // TODO: dyn allocate file list
- sParams.nDevNeighbourhood = atoi(opt.nDevNeighbourhood->answer);
-
- sParams.dumpFile = opt.dumpFile->answer;
- sParams.outputSeries = opt.outputSeries->answer;
-
- sParams.parcelSizeFile = opt.parcelSizeFile->answer;
-
- sParams.discountFactor = atof(opt.discountFactor->answer);
-
- // always 1 if not stochastic 2
- sParams.sortProbs = 1;
-
- // TODO: implement real switching of algorithm
- if (!strcmp(opt.algorithm->answer, "deterministic"))
- sParams.nAlgorithm = _N_ALGORITHM_DETERMINISTIC;
- else if (!strcmp(opt.algorithm->answer, "stochastic1"))
- sParams.nAlgorithm = _N_ALGORITHM_STOCHASTIC_I;
- else if (!strcmp(opt.algorithm->answer, "stochastic2"))
- sParams.nAlgorithm = _N_ALGORITHM_STOCHASTIC_II;
-
- if (sParams.nAlgorithm == _N_ALGORITHM_STOCHASTIC_I) {
- // TODO: add check of filled answer
- sParams.dProbWeight = atof(opt.dProbWeight->answer);
- sParams.dDevPersistence = atof(opt.dDevPersistence->answer);
-
- sParams.giveUpRatio = atof(opt.giveUpRatio->answer);
- }
- else if (sParams.nAlgorithm == _N_ALGORITHM_STOCHASTIC_II) {
-
- int parsedOK, i;
- FILE *fp;
- char inBuff[N_MAXREADINLEN];
- char *pPtr;
-
- fprintf(stdout, "reading probability lookup\n");
- sParams.probLookupFile = opt.probLookupFile->answer;
-
- fp = fopen(sParams.probLookupFile, "r");
- if (fp) {
- parsedOK = 0;
- if (fgets(inBuff, N_MAXREADINLEN, fp)) {
- if (inBuff[0] == ',') {
- sParams.nProbLookup = atoi(inBuff + 1);
- if (sParams.nProbLookup > 0) {
- sParams.adProbLookup =
- (double *)malloc(sizeof(double) *
- sParams.nProbLookup);
- if (sParams.adProbLookup) {
- parsedOK = 1;
- i = 0;
- while (parsedOK && i < sParams.nProbLookup) {
- parsedOK = 0;
- if (fgets(inBuff, N_MAXREADINLEN, fp)) {
- if (pPtr = strchr(inBuff, ',')) {
- parsedOK = 1;
- sParams.adProbLookup[i] =
- atof(pPtr + 1);
- G_debug(3,
- "probability lookup table: i=%d, i/(n-1)=%f, result=%f",
- i,
- i * 1.0 /
- (sParams.nProbLookup - 1),
- sParams.adProbLookup[i]);
- }
- }
- i++;
- }
- }
- }
- }
- }
- if (!parsedOK) {
- fprintf(stderr,
- "error parsing probLookup file '%s'...exiting\n",
- sParams.probLookupFile);
- return 0;
- }
- fclose(fp);
- }
- else {
- perror("The following error occurred");
- fprintf(stderr, "error opening probLookup file '%s'...exiting\n",
- sParams.probLookupFile);
- return 0;
- }
-
- sParams.sortProbs = atoi(opt.sortProbs->answer);
- sParams.patchFactor = atof(opt.patchFactor->answer);
- sParams.patchMean = atof(opt.patchMean->answer);
- sParams.patchRange = atof(opt.patchRange->answer);
- sParams.numNeighbors = atoi(opt.numNeighbors->answer);
- // TODO: convert to options or flag: 1: uniform distribution 2: based on dev. proba.
- sParams.seedSearch = atoi(opt.seedSearch->answer);
- sParams.devPressureApproach = atoi(opt.devPressureApproach->answer);
- if (sParams.devPressureApproach != 1) {
- sParams.alpha = atof(opt.alpha->answer);
- sParams.scalingFactor = atof(opt.scalingFactor->answer);
- }
-
- sParams.num_Regions = atoi(opt.num_Regions->answer);
- sParams.indexFile = opt.indexFile->answer;
- sParams.controlFileAll = opt.controlFileAll->answer;
- }
-
- if (sParams.num_Regions > 1)
- readDevPotParams(&sParams, opt.devpotParamsFile->answer);
-
- readDevDemand(&sParams);
- /* allocate memory */
- if (buildLandscape(&sLandscape, &sParams)) {
- /* read data */
- if (readData(&sLandscape, &sParams)) {
- readData4AdditionalVariables(&sLandscape, &sParams);
- readIndexData(&sLandscape, &sParams);
- if (readParcelSizes(&sLandscape, &sParams)) {
- //testDevPressure(&sLandscape, &sParams);
- /* do calculation and dump result */
- //one region, enable the following
- //updateMap(&sLandscape, &sParams);
-
- //multiple regions
- updateMapAll(&sLandscape, &sParams);
- }
- else {
- fprintf(stderr, "error in readParcelSizes()\n");
- }
- }
- else {
- fprintf(stderr, "error in readData()\n");
- }
- /* could put in routines to free memory, but OS will garbage collect anyway */
- }
- else {
- fprintf(stderr, "error in buildLandscape()\n");
- }
-
- return EXIT_SUCCESS;
-}
-
-int getUnDevIndex(t_Landscape * pLandscape)
-{
- float p = rand() / (double)RAND_MAX;
- int i;
-
- for (i = 0; i < pLandscape->undevSites; i++) {
- if (p < pLandscape->asUndev[i].cumulProb) {
- return i;
- }
- }
- // TODO: returning at least something but should be something more meaningful
- return 0;
-}
-
-int getUnDevIndex1(t_Landscape * pLandscape, int regionID)
-{
- float p = rand() / (double)RAND_MAX;
-
- G_verbose_message(_("getUnDevIndex1: regionID=%d, num_undevSites=%d, p=%f"),
- regionID, pLandscape->num_undevSites[regionID], p);
- int i;
-
- for (i = 0; i < pLandscape->num_undevSites[regionID]; i++) {
- if (p < pLandscape->asUndevs[regionID][i].cumulProb) {
- G_verbose_message(_("getUnDevIndex1: cumulProb=%f"),
- pLandscape->asUndevs[regionID][i].cumulProb);
- return i;
- }
- }
- // TODO: returning at least something but should be something more meaningful
- return 0;
-}
-
-void findAndSortProbsAll(t_Landscape * pLandscape, t_Params * pParams,
- int step)
-{
- int i, lookupPos;
- t_Cell *pThis;
- int id;
-
- /* update calcs */
- fprintf(stdout, "\t\trecalculating probabilities\n");
- for (i = 0; i < pParams->num_Regions; i++) {
- pLandscape->num_undevSites[i] = 0;
- }
- float val = 0.0;
-
- for (i = 0; i < pLandscape->totalCells; i++) {
- pThis = &(pLandscape->asCells[i]);
- if (pThis->nCellType == _CELL_VALID) {
- if (pThis->bUndeveloped) {
- if (pThis->consWeight > 0.0) {
- id = pThis->index_region - 1;
- if (pThis->index_region == -9999)
- continue;
- /* note that are no longer just storing the logit value,
- but instead the probability
- (allows consWeight to affect sort order) */
-
- if (id < 0 || id >= pParams->num_Regions)
- G_fatal_error(_("Index of region %d is out of range [0, %d] (region is %d)"),
- id, pParams->num_Regions - 1,
- pThis->index_region);
-
- if (pLandscape->num_undevSites[id] >=
- pLandscape->asUndevs[id].size())
- pLandscape->asUndevs[id].resize(pLandscape->
- asUndevs[id].size() *
- 2);
-
- pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].
- cellID = i;
- val = getDevProbability(pThis, pParams);
- pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].
- logitVal = val;
- G_verbose_message("logit value %f", val);
- pThis->devProba = val;
- /* lookup table of probabilities is applied before consWeight */
- if (pParams->nAlgorithm == _N_ALGORITHM_STOCHASTIC_II) {
- /* replace with value from lookup table */
- lookupPos = (int)(pLandscape->asUndevs[id]
- [pLandscape->num_undevSites[id]].
- logitVal * (pParams->nProbLookup -
- 1));
- if (lookupPos >= pParams->nProbLookup ||
- lookupPos < 0)
- G_fatal_error
- ("lookup position (%d) out of range [0, %d]",
- lookupPos, pParams->nProbLookup - 1);
- pLandscape->asUndevs[id][pLandscape->
- num_undevSites[id]].
- logitVal = pParams->adProbLookup[lookupPos];
- }
- pThis->devProba =
- pLandscape->asUndevs[id][pLandscape->
- num_undevSites[id]].logitVal;
- // discount by a conservation factor
- pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].logitVal *= pThis->consWeight;
- /* need to store this to put correct elements near top of list */
- pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].bUntouched = pThis->bUntouched;
- if (pLandscape->asUndevs[id]
- [pLandscape->num_undevSites[id]].logitVal > 0.0) {
- /* only add one more to the list if have a positive probability */
- pLandscape->num_undevSites[id]++;
- }
- }
- }
- }
- }
- /* downweight the devPressure if necessary (do not do in first step) */
- /* doing it here means that last time step values have full weight */
- if (pParams->nAlgorithm == _N_ALGORITHM_STOCHASTIC_I) {
- if (pParams->dDevPersistence < 1.0) {
- fprintf(stdout, "\t\tdownweighting development pressure\n");
-
- for (i = 0; i < pLandscape->totalCells; i++) {
- pThis = &(pLandscape->asCells[i]);
- if (pThis->nCellType == _CELL_VALID) {
- /* only need to bother downweighting on cells that can still convert */
- if (pThis->bUndeveloped) {
- pThis->devPressure =
- (int)((double)pThis->devPressure *
- pParams->dDevPersistence);
- }
- }
- }
- }
- }
- /* sort */
- /* can only be zero for algorithm=stochastic_ii */
- if (pParams->sortProbs) {
- fprintf(stdout, "\t\tsorting %d unconserved undeveloped sites\n",
- pLandscape->undevSites);
- qsort(pLandscape->asUndev, pLandscape->undevSites, sizeof(t_Undev),
- undevCmpReverse);
- }
- else {
- fprintf(stdout, "\t\tskipping sort as choosing cells randomly\n");
- }
- //calculate cumulative probability // From Wenwu Tang
- int j;
-
- for (j = 0; j < pParams->num_Regions; j++) {
- double sum = pLandscape->asUndevs[j][0].logitVal;
-
- for (i = 1; i < pLandscape->num_undevSites[j]; i++) {
- pLandscape->asUndevs[j][i].cumulProb =
- pLandscape->asUndevs[j][i - 1].cumulProb +
- pLandscape->asUndevs[j][i].logitVal;
- sum = sum + pLandscape->asUndevs[j][i].logitVal;
- }
- for (i = 0; i < pLandscape->num_undevSites[j]; i++) {
- pLandscape->asUndevs[j][i].cumulProb =
- pLandscape->asUndevs[j][i].cumulProb / sum;
- }
- }
-}
Added: grass-addons/grass7/raster/r.futures/r.futures.calib/Makefile
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.calib/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.calib/Makefile 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.futures.calib
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.html 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,113 @@
+<h2>DESCRIPTION</h2>
+Module <em>r.futures.calibration</em> is part of <a href="r.futures.html">FUTURES</a>,
+land change model.
+It is used for calibrating certain input variables for patch growing
+algorithm <em>r.futures.pga</em>, specifically patch size and compactness parameters.
+The calibration process is conducted to match observed urban growth patterns
+to those simulated by the model, including the sizes and shapes of new development.
+The calibration is achieved by varying the values of the patch parameters,
+comparing the distribution of simulated patch sizes to those observed
+for the reference period, and choosing the values that provide the closest match.
+For the details about calibration see below.
+
+
+<h3>Patch size</h3>
+As part of the calibration process, module <em>r.futures.calibration</em>
+produces patch size distribution file specified in <b>patch_sizes</b> parameter,
+which contains sizes (in cells) of all new patches observed
+in the reference period. The format of this file is one patch size per line.
+FUTURES uses this file to determine the size of the simulated patches.
+Often the length of the reference time period does not match
+the time period which we are trying to simulate.
+We use the <b>discount factor</b> to alter the size of simulated patches
+so that after the reference period they closely match the observed patterns.
+During the simulation, this factor is multiplied by the patch sizes listed in the patch size file.
+The values of <b>discount factor</b> can vary between 0 and 1,
+for example value 0.6 was used by Meentemeyer et al. 2013.
+
+<h3>Patch compactness</h3>
+The shapes of patches simulated by FUTURES are governed
+by the patch compactness parameter (Meentemeyer et al. 2013, Eq. 1).
+This variable doesn't represent actual patch compactness, it is rather
+an adjustable scaling factor that controls patch compactness
+through a distance decay effect.
+By specifying the mean and range of this parameter in module
+<em>r.futures.pga</em>, we allow for variation in patch shape.
+As the value of the parameter increases, patches become more compact.
+Calibration is achieved by varying the values specified in <b>compactness_mean</b>
+and <b>compactness_range</b> and comparing the distribution
+of the simulated patch compactness (computed as patch perimeter / (2 * sqrt(pi * area)))
+to those observed for the reference period.
+Meentemeyer et al. 2013 used mean 0.4 and range 0.08.
+
+<h3>Calibration input and output</h3>
+Calibration requires the development binary raster in the beginning
+and end of the reference period (<b>development_start</b> and <b>development_end</b>)
+to derive the patch sizes and compactness.
+It is possible to set the minimum number of cells of a patch in <b>patch_threshold</b>
+to ignore too small patches.
+For each combination of values provided in <b>compactness_mean</b>,
+<b>compactness_range</b> and <b>discount_factor</b>, it runs
+module <em>r.futures.pga</em> which creates new development pattern.
+From this new simulated development, patch characteristics are derived
+and compared with the observed characteristics by histogram comparison
+and an error (histogram distance) is computed.
+Since <em>r.futures.pga</em> is a stochastic module, multiple
+runs (specified in <b>repeat</b>) are recommended, the error is then averaged.
+Calibration results are saved in a csv file specified in <b>calibration_results</b>:
+
+<pre>
+input_discount_factor area_distance input_compactness_mean input_compactness_range compactness_distance
+0.1 1.01541178435 0.1 0.02 3.00000005937
+0.2 1.26578803108 0.1 0.02 4.12442780529
+0.3 1.17631210026 0.1 0.02 3.86904462396
+0.4 2.31700278644 0.1 0.02 15.0569602795
+0.5 1.08655152036 0.1 0.02 3.72484862687
+0.6 2.97628078734 0.1 0.02 21.6358616001
+0.7 3.61632549044 0.1 0.02 25.4492265706
+0.8 2.72789958233 0.1 0.02 18.1083820007
+0.9 2.45915297845 0.1 0.02 18.4500322711
+0.1 1.05473877995 0.1 0.04 3.09321560218
+...
+</pre>
+
+Optimal values can be found by visual examination of the second and fifth columns.
+<p>
+Providing too many values in <b>compactness_mean</b>,
+<b>compactness_range</b> and <b>discount_factor</b> results in very long computation.
+Therefore it is recommended to run <em>r.futures.calibration</em>
+on high-end computers, with more processes running in parallel using <b>nprocs</b> parameter.
+Also, it can be run on smaller region, under the assumption
+that patch sizes and shapes are close to being consistent across the entire study area.
+
+<p>
+For all other parameters not mentioned above, please refer to
+<em>r.futures.pga</em> documentation.
+
+<h2>EXAMPLES</h2>
+
+<h2>SEE ALSO</h2>
+
+<a href="r.futures.html">FUTURES</a>,
+<em><a href="r.futures.pga.html">r.futures.pga</a></em>,
+<em><a href="r.futures.devpressure.html">r.futures.devpressure</a></em>,
+<em><a href="r.sample.category.html">r.sample.category</a></em>
+
+<h2>REFERENCES</h2>
+<ul>
+<li>
+ Meentemeyer, R. K., Tang, W., Dorning, M. A., Vogler, J. B., Cunniffe, N. J., & Shoemaker, D. A. (2013).
+ <a href="http://dx.doi.org/10.1080/00045608.2012.707591">FUTURES: Multilevel Simulations of Emerging
+ Urban-Rural Landscape Structure Using a Stochastic Patch-Growing Algorithm</a>.
+ Annals of the Association of American Geographers, 103(4), 785-807.
+<li>Dorning, M. A., Koch, J., Shoemaker, D. A., & Meentemeyer, R. K. (2015).
+ <a href="http://dx.doi.org/10.1016/j.landurbplan.2014.11.011">Simulating urbanization scenarios reveals
+ tradeoffs between conservation planning strategies</a>.
+ Landscape and Urban Planning, 136, 28-39.</li>
+</ul>
+
+<h2>AUTHOR</h2>
+
+Anna Petrasova, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a><br>
+
+<p><i>Last changed: $Date: 2015-02-12 14:32:45 -0500 (Thu, 12 Feb 2015) $</i>
Added: grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.py
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.py (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.calib/r.futures.calib.py 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,453 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+##############################################################################
+#
+# MODULE: r.futures.calib
+#
+# AUTHOR(S): Anna Petrasova (kratochanna gmail.com)
+#
+# PURPOSE: FUTURES patches calibration tool
+#
+# COPYRIGHT: (C) 2015 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+##############################################################################
+
+#%module
+#% description: Script for calibrating patch characteristics used as input to r.futures.pga
+#% keyword: raster
+#% keyword: patch
+#%end
+#%option G_OPT_R_INPUT
+#% key: development_start
+#% label: Name of input binary raster map representing development in the beginning
+#% description: Raster map of developed areas (=1), undeveloped (=0) and excluded (no data)
+#% guisection: Calibration
+#%end
+#%option G_OPT_R_INPUT
+#% key: development_end
+#% label: Name of input binary raster map representing development in the end
+#% description: Raster map of developed areas (=1), undeveloped (=0) and excluded (no data)
+#% guisection: Calibration
+#%end
+#%option
+#% type: integer
+#% key: repeat
+#% description: How many times is the simulation repeated
+#% required: yes
+#% answer: 10
+#% guisection: Calibration
+#%end
+#%option
+#% key: compactness_mean
+#% type: double
+#% description: Patch compactness mean to be tested
+#% required: yes
+#% multiple: yes
+#% guisection: Calibration
+#%end
+#%option
+#% type: double
+#% key: compactness_range
+#% description: Patch compactness range to be tested
+#% required: yes
+#% multiple: yes
+#% guisection: Calibration
+#%end
+#%option
+#% type: double
+#% key: discount_factor
+#% description: Patch size discount factor
+#% key: discount_factor
+#% multiple: yes
+#% required: yes
+#% guisection: Calibration
+#%end
+#%option
+#% type: double
+#% key: patch_threshold
+#% description: Minimum size of a patch in meters squared
+#% required: yes
+#% answer: 0
+#% guisection: Calibration
+#%end
+#%option G_OPT_F_OUTPUT
+#% key: patch_sizes
+#% description: Output file with patch sizes
+#% required: yes
+#% guisection: Calibration
+#%end
+#%option G_OPT_F_OUTPUT
+#% key: calibration_results
+#% description: Output file with calibration results
+#% required: yes
+#% guisection: Calibration
+#%end
+#%option
+#% key: nprocs
+#% type: integer
+#% description: Number of parallel processes
+#% required: yes
+#% answer: 1
+#% guisection: Calibration
+#%end
+#%option G_OPT_R_INPUT
+#% key: development_pressure
+#% required: yes
+#% description: Files containing the information to read in
+#% guisection: PGA
+#%end
+#%option G_OPT_R_INPUT
+#% key: cons_weight
+#% required: no
+#% label: Name of raster map representing development potential constraint weight for scenarios
+#% description: Values must be between 0 and 1, 1 means no constraint
+#% guisection: PGA
+#%end
+#%option G_OPT_R_INPUTS
+#% key: predictors
+#% required: yes
+#% multiple: yes
+#% description: Names of predictor variable raster maps
+#% guisection: PGA
+#%end
+#%option
+#% key: n_dev_neighbourhood
+#% type: integer
+#% description: Size of square used to recalculate development pressure
+#% required: yes
+#% guisection: PGA
+#%end
+#%option G_OPT_F_INPUT
+#% key: devpot_params
+#% required: yes
+#% multiple: yes
+#% label: Development potential parameters for each region
+#% description: Each line should contain region ID followed by parameters. Values are separated by whitespace (spaces or tabs). First line is ignored, so it can be used for header
+#% guisection: PGA
+#%end
+#%option G_OPT_F_INPUT
+#% key: incentive_table
+#% required: yes
+#% description: File containing incentive lookup table (infill vs. sprawl)
+#% guisection: PGA
+#%end
+#%option
+#% key: num_neighbors
+#% type: integer
+#% required: yes
+#% multiple: no
+#% options: 4,8
+#% description: The number of neighbors to be used for patch generation (4 or 8)
+#% guisection: PGA
+#%end
+#%option
+#% key: seed_search
+#% type: integer
+#% required: yes
+#% multiple: no
+#% options: 1,2
+#% description: The way that the location of a seed is determined
+#% guisection: PGA
+#%end
+#%option
+#% key: development_pressure_approach
+#% type: string
+#% required: yes
+#% multiple: no
+#% options: occurrence,gravity,kernel
+#% answer: gravity
+#% description: Approaches to derive development pressure
+#% guisection: PGA
+#%end
+#%option
+#% key: gamma
+#% type: double
+#% required: yes
+#% multiple: no
+#% description: Influence of distance between neighboring cells
+#% guisection: PGA
+#%end
+#%option
+#% key: scaling_factor
+#% type: double
+#% required: yes
+#% multiple: no
+#% description: Scaling factor
+#% guisection: PGA
+#%end
+#%option
+#% key: num_regions
+#% type: integer
+#% required: yes
+#% multiple: no
+#% description: Number of sub-regions (e.g., counties) to be simulated
+#% guisection: PGA
+#%end
+#%option G_OPT_R_INPUT
+#% key: subregions
+#% required: yes
+#% description: Raster map of subregions with categories starting with 1
+#% guisection: PGA
+#%end
+#%option G_OPT_F_INPUT
+#% key: demand
+#% required: yes
+#% description: Control file with number of cells to convert
+#% guisection: PGA
+#%end
+
+
+import sys
+import os
+import atexit
+import numpy as np
+import tempfile
+from multiprocessing import Process, Queue
+
+import grass.script.core as gcore
+import grass.script.raster as grast
+import grass.script.utils as gutils
+
+
+TMP = []
+TMPFILE = None
+
+
+def cleanup(tmp=None):
+ if tmp:
+ maps = tmp
+ else:
+ maps = TMP
+ gcore.run_command('g.remove', flags='f', type=['raster', 'vector'], name=maps)
+ gutils.try_remove(TMPFILE)
+
+
+def run_one_combination(repeat, development_start, compactness_mean, compactness_range,
+ discount_factor, patches_file, fut_options, threshold,
+ hist_bins_area_orig, hist_range_area_orig, hist_bins_compactness_orig,
+ hist_range_compactness_orig, cell_size, histogram_area_orig, histogram_compactness_orig,
+ tmp_name, queue):
+ TMP_PROCESS = []
+ # unique name, must be sql compliant
+ suffix = (str(discount_factor) + str(compactness_mean) + str(compactness_range)).replace('.', '')
+ simulation_dev_end = tmp_name + 'simulation_dev_end_' + suffix
+ simulation_dev_diff = tmp_name + 'simulation_dev_diff' + suffix
+ tmp_patch_vect = tmp_name + 'tmp_patch_vect' + suffix
+ tmp_patch_vect2 = tmp_name + 'tmp_patch_vect2' + suffix
+ TMP_PROCESS.append(simulation_dev_diff)
+ TMP_PROCESS.append(simulation_dev_diff)
+ TMP_PROCESS.append(tmp_patch_vect)
+ TMP_PROCESS.append(tmp_patch_vect2)
+
+ sum_dist_area = 0
+ sum_dist_compactness = 0
+ for i in range(repeat):
+ gcore.message(_("Running FUTURES simulation {i}/{repeat}...".format(i=i + 1, repeat=repeat)))
+ run_simulation(development_start=development_start, development_end=simulation_dev_end,
+ compactness_mean=compactness_mean, compactness_range=compactness_range,
+ discount_factor=discount_factor, patches_file=patches_file, fut_options=fut_options)
+ new_development(simulation_dev_end, simulation_dev_diff)
+
+ temp_file = tempfile.NamedTemporaryFile(delete=False)
+ temp_file.close()
+ patch_analysis(simulation_dev_diff, threshold, tmp_patch_vect, tmp_patch_vect2, temp_file.name)
+ sim_hist_area, sim_hist_compactness = create_histograms(temp_file.name, hist_bins_area_orig, hist_range_area_orig,
+ hist_bins_compactness_orig, hist_range_compactness_orig, cell_size)
+ os.remove(temp_file.name)
+
+ sum_dist_area += compare_histograms(histogram_area_orig, sim_hist_area)
+ sum_dist_compactness += compare_histograms(histogram_compactness_orig, sim_hist_compactness)
+
+ mean_dist_area = sum_dist_area / repeat
+ mean_dist_compactness = sum_dist_compactness / repeat
+
+ data = {}
+ data['input_discount_factor'] = discount_factor
+ data['input_compactness_mean'] = compactness_mean
+ data['input_compactness_range'] = compactness_range
+ data['area_distance'] = mean_dist_area
+ data['compactness_distance'] = mean_dist_compactness
+ queue.put(data)
+ cleanup(tmp=TMP_PROCESS)
+
+
+def run_simulation(development_start, development_end, compactness_mean, compactness_range, discount_factor, patches_file, fut_options):
+ parameters = dict(patch_mean=compactness_mean, patch_range=compactness_range,
+ discount_factor=discount_factor, patch_sizes=patches_file,
+ developed=development_start)
+ futures_parameters = dict(development_pressure=fut_options['development_pressure'],
+ predictors=fut_options['predictors'], n_dev_neighbourhood=fut_options['n_dev_neighbourhood'],
+ devpot_params=fut_options['devpot_params'], incentive_table=fut_options['incentive_table'],
+ num_neighbors=fut_options['num_neighbors'], seed_search=fut_options['seed_search'],
+ development_pressure_approach=fut_options['development_pressure_approach'], gamma=fut_options['gamma'],
+ scaling_factor=fut_options['scaling_factor'], num_regions=fut_options['num_regions'],
+ subregions=fut_options['subregions'], demand=fut_options['demand'],
+ output=development_end)
+ parameters.update(futures_parameters)
+ for not_required in ('cons_weight',):
+ if options[not_required]:
+ parameters.update({not_required: options[not_required]})
+
+ gcore.run_command('r.futures.pga', flags='s', overwrite=True, **parameters)
+
+
+def diff_development(development_start, development_end, subregions, development_diff):
+ grast.mapcalc(exp="{res} = if({subregions} && {dev_end} && (isnull({dev_start}) ||| !{dev_start}), 1, null())".format(res=development_diff, subregions=subregions,
+ dev_end=development_end, dev_start=development_start), overwrite=True, quiet=True)
+
+
+def new_development(development_end, development_diff):
+ grast.mapcalc(exp="{res} = if({dev_end} > 0, 1, null())".format(res=development_diff,
+ dev_end=development_end), overwrite=True, quiet=True)
+
+
+def patch_analysis(development_diff, threshold, tmp_vector_patches, tmp_vector_patches2, output_file):
+ gcore.run_command('r.to.vect', input=development_diff, output=tmp_vector_patches2, type='area', overwrite=True, quiet=True)
+ gcore.run_command('v.clean', input=tmp_vector_patches2, output=tmp_vector_patches, tool='rmarea', threshold=threshold, quiet=True, overwrite=True)
+ gcore.run_command('v.db.addcolumn', map=tmp_vector_patches, columns="area double precision,perimeter double precision", quiet=True)
+ gcore.run_command('v.to.db', map=tmp_vector_patches, option='area', column='area', units='meters', quiet=True)
+ gcore.run_command('v.to.db', map=tmp_vector_patches, option='perimeter', column='perimeter', units='meters', quiet=True)
+ gcore.run_command('v.db.select', map=tmp_vector_patches, columns=['area', 'perimeter'],
+ flags='c', separator='space', file_=output_file, overwrite=True, quiet=True)
+
+
+def create_histograms(input_file, hist_bins_area_orig, hist_range_area_orig, hist_bins_compactness_orig, hist_range_compactness_orig, cell_size):
+ area, perimeter = np.loadtxt(fname=input_file, unpack=True)
+ compact = compactness(area, perimeter)
+ histogram_area, _edges = np.histogram(area / cell_size, bins=hist_bins_area_orig,
+ range=hist_range_area_orig, density=True)
+ histogram_area = histogram_area * 100
+ histogram_compactness, _edges = np.histogram(compact, bins=hist_bins_compactness_orig,
+ range=hist_range_compactness_orig, density=True)
+ histogram_compactness = histogram_compactness * 100
+ return histogram_area, histogram_compactness
+
+
+def write_patches_file(vector_patches, cell_size, output_file):
+ gcore.run_command('v.db.select', map=vector_patches, columns='area',
+ flags='c', separator='space', file_=output_file, quiet=True)
+ areas = np.loadtxt(fname=output_file)
+ areas = np.round(areas / cell_size)
+ np.savetxt(fname=output_file, X=areas.astype(int), fmt='%u')
+
+
+def compare_histograms(hist1, hist2):
+ """
+ >>> hist1, edg = np.histogram(np.array([1, 1, 2, 2.5, 2.4]), bins=3, range=(0, 6))
+ >>> hist2, edg = np.histogram(np.array([1, 1, 3 ]), bins=3, range=(0, 6))
+ >>> compare_histograms(hist1, hist2)
+ 0.5
+ """
+ mask = np.logical_not(np.logical_or(hist1, hist2))
+ hist1 = np.ma.masked_array(hist1, mask=mask)
+ hist2 = np.ma.masked_array(hist2, mask=mask)
+ res = 0.5 * np.sum(np.power(hist1 - hist2, 2) / (hist1.astype(float) + hist2))
+ return res
+
+
+def compactness(area, perimeter):
+ return perimeter / (2 * np.sqrt(np.pi * area))
+
+
+def main():
+ dev_start = options['development_start']
+ dev_end = options['development_end']
+ repeat = int(options['repeat'])
+ compactness_means = [float(each) for each in options['compactness_mean'].split(',')]
+ compactness_ranges = [float(each) for each in options['compactness_range'].split(',')]
+ discount_factors = [float(each) for each in options['discount_factor'].split(',')]
+ patches_file = options['patch_sizes']
+ threshold = float(options['patch_threshold'])
+ # v.clean removes size <= threshold, we want to keep size == threshold
+ threshold -= 1e-6
+
+ # compute cell size
+ region = gcore.region()
+ res = (region['nsres'] + region['ewres'])/2.
+ coeff = float(gcore.parse_command('g.proj', flags='g')['meters'])
+ cell_size = res * res * coeff * coeff
+
+ tmp_name = 'tmp_futures_calib_' + str(os.getpid()) + '_'
+ global TMP, TMPFILE
+
+ orig_patch_diff = tmp_name + 'orig_patch_diff'
+ TMP.append(orig_patch_diff)
+ tmp_patch_vect = tmp_name + 'tmp_patch_vect'
+ tmp_patch_vect2 = tmp_name + 'tmp_patch_vect2'
+ TMP.append(tmp_patch_vect)
+ TMP.append(tmp_patch_vect2)
+ temp_file = tempfile.NamedTemporaryFile(delete=False)
+ temp_file.close()
+ TMPFILE = temp_file.name
+
+ gcore.message(_("Analyzing original patches..."))
+ diff_development(dev_start, dev_end, options['subregions'], orig_patch_diff)
+ patch_analysis(orig_patch_diff, threshold, tmp_patch_vect, tmp_patch_vect2, temp_file.name)
+ area, perimeter = np.loadtxt(fname=temp_file.name, unpack=True)
+ compact = compactness(area, perimeter)
+ write_patches_file(tmp_patch_vect, cell_size, patches_file)
+
+ # area histogram
+ area = area / cell_size
+ bin_width = 1. # automatic ways to determine bin width do not perform well in this case
+ hist_bins_area_orig = int(np.ptp(area) / bin_width)
+ hist_range_area_orig = (np.min(area), np.max(area))
+ histogram_area_orig, _edges = np.histogram(area, bins=hist_bins_area_orig,
+ range=hist_range_area_orig, density=True)
+ histogram_area_orig = histogram_area_orig * 100 # to get percentage for readability
+
+ # compactness histogram
+ bin_width = 0.1
+ hist_bins_compactness_orig = int(np.ptp(compact) / bin_width)
+ hist_range_compactness_orig = (np.min(compact), np.max(compact))
+ histogram_compactness_orig, _edges = np.histogram(compact, bins=hist_bins_compactness_orig,
+ range=hist_range_compactness_orig, density=True)
+ histogram_compactness_orig = histogram_compactness_orig * 100 # to get percentage for readability
+
+ nprocs = int(options['nprocs'])
+ count = 0
+ proc_count = 0
+ queue_list = []
+ proc_list = []
+ num_all = len(compactness_means) * len(compactness_ranges) * len(discount_factors)
+ with open(options['calibration_results'], 'a') as f:
+ f.write(' '.join(['input_discount_factor', 'area_distance',
+ 'input_compactness_mean', 'input_compactness_range',
+ 'compactness_distance']))
+ f.write('\n')
+ for com_mean in compactness_means:
+ for com_range in compactness_ranges:
+ for discount_factor in discount_factors:
+ count += 1
+ q = Queue()
+ p = Process(target=run_one_combination,
+ args=(repeat, dev_start, com_mean, com_range,
+ discount_factor, patches_file, options, threshold,
+ hist_bins_area_orig, hist_range_area_orig, hist_bins_compactness_orig,
+ hist_range_compactness_orig, cell_size, histogram_area_orig, histogram_compactness_orig,
+ tmp_name, q))
+ p.start()
+ queue_list.append(q)
+ proc_list.append(p)
+ proc_count += 1
+ if proc_count == nprocs or count == num_all:
+ for i in range(proc_count):
+ proc_list[i].join()
+ data = queue_list[i].get()
+ f.write(' '.join([str(data['input_discount_factor']), str(data['area_distance']),
+ str(data['input_compactness_mean']), str(data['input_compactness_range']),
+ str(data['compactness_distance'])]))
+ f.write('\n')
+ proc_count = 0
+ proc_list = []
+ queue_list = []
+
+if __name__ == "__main__":
+ options, flags = gcore.parser()
+ atexit.register(cleanup)
+ sys.exit(main())
Added: grass-addons/grass7/raster/r.futures/r.futures.devpressure/Makefile
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.devpressure/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.devpressure/Makefile 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,7 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.futures.devpressure
+
+include $(MODULE_TOPDIR)/include/Make/Script.make
+
+default: script
Added: grass-addons/grass7/raster/r.futures/r.futures.devpressure/devpressure_example.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.futures/r.futures.devpressure/devpressure_example.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Added: grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.html 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,81 @@
+<h2>DESCRIPTION</h2>
+Module <em>r.futures.devpressure</em> is part of <a href="r.futures.html">FUTURES</a>,
+land change model.
+It computes development pressure, a predictor based on number of neighboring
+developed cells within search distance, weighted by distance.
+The development pressure variable plays a special role in the model,
+allowing for a feedback between predicted change and change in subsequent steps.
+
+<p>
+There are multiple options for calculating development pressure in the model.
+<ul>
+<li><b>occurrence</b>: simple count of the number of developed
+cells within the specified window size</li>
+<li><b>gravity</b>: defined as scaling_factor / pow(distance, gamma)</li>
+<li><b>kernel</b>: defined as scaling_factor * exp(-2 * distance / gamma)</li>
+</ul>
+<p>
+Coefficient <b>gamma</b> is the coefficient controlling the influence of distance.
+The best value for gamma is chosen during the model selection process (TODO).
+
+<p>
+The input raster is binary, value 1 represents development, value 0 no development.
+The neighborhood <b>size</b> describes the size of the moving
+window in cells (2 * size + 1).
+
+<p>
+<center>
+<img src="devpressure_example.png">
+<p>
+Figure: Effect of parameters <b>size</b> and <b>gamma</b>:
+a) initial development,
+b) gamma = 2, size = 10,
+c) gamma = 1.5, size = 10,
+d) gamma = 1.5, size = 20.
+</center>
+
+
+
+<h2>NOTES</h2>
+By default NULL values are propagated.
+This means that the edges of resulting map will be NULLs. To avoid that,
+flag <b>n</b> internally enlarges the computational region,
+converts NULLs to zeros, performs the computation
+and then patches back the original NULL values.
+
+<p>
+Module <em>r.futures.devpressure</em>, although written for FUTURES model,
+is general enough to be used for different applications where distance pressure
+ can be described with the functions above.
+
+
+<h2>EXAMPLES</h2>
+<div class="code"><pre>
+r.futures.devpressure input=developed output=pressure gamma=1.5 size=10 method=gravity
+</pre></div>
+
+<h2>SEE ALSO</h2>
+
+<a href="r.futures.html">FUTURES</a>,
+<em><a href="r.futures.pga.html">r.futures.pga</a></em>,
+<em><a href="r.futures.calib.html">r.futures.calib</a></em>,
+<em><a href="r.sample.category.html">r.sample.category</a></em>
+
+<h2>REFERENCES</h2>
+<ul>
+<li>
+ Meentemeyer, R. K., Tang, W., Dorning, M. A., Vogler, J. B., Cunniffe, N. J., & Shoemaker, D. A. (2013).
+ <a href="http://dx.doi.org/10.1080/00045608.2012.707591">FUTURES: Multilevel Simulations of Emerging
+ Urban-Rural Landscape Structure Using a Stochastic Patch-Growing Algorithm</a>.
+ Annals of the Association of American Geographers, 103(4), 785-807.
+<li>Dorning, M. A., Koch, J., Shoemaker, D. A., & Meentemeyer, R. K. (2015).
+ <a href="http://dx.doi.org/10.1016/j.landurbplan.2014.11.011">Simulating urbanization scenarios reveals
+ tradeoffs between conservation planning strategies</a>.
+ Landscape and Urban Planning, 136, 28-39.</li>
+</ul>
+
+<h2>AUTHOR</h2>
+
+Anna Petrasova, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a><br>
+
+<p><i>Last changed: $Date: 2015-02-12 14:32:45 -0500 (Thu, 12 Feb 2015) $</i>
Added: grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.py
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.py (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.devpressure/r.futures.devpressure.py 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,187 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+#
+##############################################################################
+#
+# MODULE: r.futures.devpressure
+#
+# AUTHOR(S): Anna Petrasova (kratochanna gmail.com)
+#
+# PURPOSE: FUTURES development pressure computation
+#
+# COPYRIGHT: (C) 2015 by the GRASS Development Team
+#
+# This program is free software under the GNU General Public
+# License (>=v2). Read the file COPYING that comes with GRASS
+# for details.
+#
+##############################################################################
+
+#%module
+#% description: Script for computing development pressure
+#% keyword: raster
+#% keyword: filter
+#% keyword: statistics
+#%end
+#%option G_OPT_R_INPUT
+#% description: Name of input binary raster map representing development
+#%end
+#%option G_OPT_R_OUTPUT
+#% description: Name of the output development pressure raster
+#%end
+#%option
+#% key: method
+#% type: string
+#% description: Method for computing development pressure
+#% required: yes
+#% answer: gravity
+#% options: occurrence,gravity,kernel
+#% descriptions: occurrence;number of developed cells in window;gravity;scaling_factor / pow(distance, gamma);kernel;scaling_factor * exp (-2*distance / gamma)
+#%end
+#%option
+#% key: size
+#% type: integer
+#% description: Half of neighborhood size
+#% required: yes
+#% answer: 8
+#%end
+#%option
+#% key: gamma
+#% type: double
+#% description: Coefficient controlling the influence of distance, needed for method gravity and kernel
+#% required: no
+#% answer: 1.5
+#%end
+#%option
+#% key: scaling_factor
+#% type: double
+#% description: Scaling factor needed for method gravity and kernel
+#% required: no
+#% answer: 1
+#%end
+#%flag
+#% key: n
+#% description: Do not propagate nulls
+#%end
+
+
+import os
+import sys
+import atexit
+import numpy as np
+from math import sqrt
+
+#from grass.exceptions import CalledModuleError
+import grass.script.core as gcore
+import grass.script.utils as gutils
+import grass.script.raster as grast
+
+TMPFILE = None
+TMP = []
+
+
+def cleanup():
+ if TMP:
+ gcore.run_command('g.remove', flags='f', type=['raster'], name=TMP)
+ gutils.try_remove(TMPFILE)
+
+
+def main():
+ size = int(options['size'])
+ gamma = scale = None
+ if options['gamma']:
+ gamma = float(options['gamma'])
+ if options['scaling_factor']:
+ scale = float(options['scaling_factor'])
+ input_dev = options['input']
+ output = options['output']
+ method = options['method']
+
+ if method in ('gravity', 'kernel') and (gamma is None or scale is None):
+ gcore.fatal(_("Methods gravity and kernel require options scaling_factor and gamma"))
+
+ temp_map = 'tmp_futures_devPressure_' + str(os.getpid()) + '_copy'
+ temp_map_out = 'tmp_futures_devPressure_' + str(os.getpid()) + '_out'
+ temp_map_nulls = 'tmp_futures_devPressure_' + str(os.getpid()) + '_nulls'
+ global TMP, TMPFILE
+ if flags['n']:
+ gcore.message(_("Preparing data..."))
+ region = gcore.region()
+ gcore.use_temp_region()
+ gcore.run_command('g.region', n=region['n'] + size * region['nsres'],
+ s=region['s'] - size * region['nsres'],
+ e=region['e'] + size * region['ewres'],
+ w=region['w'] - size * region['ewres'])
+ TMP.append(temp_map)
+ TMP.append(temp_map_nulls)
+ TMP.append(temp_map_out)
+ exp = "{temp_map_nulls} = if(isnull({inp}), 1, null())".format(temp_map_nulls=temp_map_nulls, inp=input_dev)
+ grast.mapcalc(exp=exp)
+ grast.mapcalc(exp="{temp} = if(isnull({inp}), 0, {inp})".format(temp=temp_map, inp=input_dev))
+ rmfilter_inp = temp_map
+ rmfilter_out = temp_map_out
+ else:
+ rmfilter_inp = input_dev
+ rmfilter_out = output
+
+ matrix = distance_matrix(size)
+ if method == 'occurrence':
+ matrix[matrix > 0] = 1
+ elif method == 'gravity':
+ with np.errstate(divide='ignore'):
+ denom = np.power(matrix, gamma)
+ matrix = scale / denom
+ matrix[denom == 0] = 0
+ else:
+ matrix_ = scale * np.exp(-2 * matrix / gamma)
+ matrix = np.where(matrix > 0, matrix_, 0)
+
+ path = gcore.tempfile()
+ global TMPFILE
+ TMPFILE = path
+
+ with open(path, 'w') as f:
+ f.write(write_filter(matrix))
+ gcore.message(_("Running development pressure filter..."))
+ gcore.run_command('r.mfilter', input=rmfilter_inp, output=rmfilter_out, filter=path)
+
+ if flags['n']:
+ gcore.run_command('g.region', n=region['n'], s=region['s'], e=region['e'], w=region['w'],)
+ grast.mapcalc(exp="{out} = if(isnull({temp_null}), {rmfilter_out}, null())".format(temp_null=temp_map_nulls,
+ rmfilter_out=rmfilter_out, out=output))
+ gcore.del_temp_region()
+
+ grast.raster_history(output)
+
+
+def distance_matrix(size):
+ matrix_size = 2 * size + 1
+ matrix = np.zeros((matrix_size, matrix_size))
+ center = size
+ for i in range(matrix_size):
+ for j in range(matrix_size):
+ dist = sqrt((i - center) * (i - center) + (j - center) * (j - center))
+ if dist <= size:
+ matrix[i, j] = dist
+ return matrix
+
+
+def write_filter(matrix):
+ filter_text = ['TITLE development pressure']
+ filter_text.append('MATRIX %s' % matrix.shape[0])
+ for i in range(matrix.shape[0]):
+ line = ''
+ for j in range(matrix.shape[0]):
+ line += str(matrix[i, j])
+ line += ' '
+ filter_text.append(line)
+ filter_text.append('DIVISOR 1')
+ filter_text.append('TYPE P')
+
+ return '\n'.join(filter_text)
+
+
+if __name__ == '__main__':
+ options, flags = gcore.parser()
+ atexit.register(cleanup)
+ sys.exit(main())
Deleted: grass-addons/grass7/raster/r.futures/r.futures.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.html 2015-07-06 18:23:20 UTC (rev 65540)
+++ grass-addons/grass7/raster/r.futures/r.futures.html 2015-07-06 18:56:57 UTC (rev 65541)
@@ -1,108 +0,0 @@
-<h2>DESCRIPTION</h2>
-
-<em>r.futures</em> is an implementation of
-FUTure Urban-Regional Environment Simulation (FUTURES)
-which is a model for multilevel simulations of emerging urban-rural
-landscape structure.
-This module uses stochastic Patch-Growing Algorithm (PGA)
-and a combination of field-based and object-based representations
-to simulate land changes.
-
-The FUTure Urban-Regional Environment Simulation (FUTURES) produces
-regional projections of landscape patterns using coupled submodels
-that integrate nonstationary drivers of land change: per capita demand,
-site suitability, and the spatial structure of conversion events.
-
-The PGA is a stochastic simulation that maps the
-allocation and spatial structure of land change using
-iterative site selection and a contextually aware region
-growing mechanism that changes cells from "undeveloped"
-to terminal "developed" states.
-
-<p>
-<center>
-<img src="r_futures.png">
-<p>
-Figure: Output map of developed areas
-</center>
-
-<p>
-<center>
-<img src="r_futures_detail.png">
-<p>
-Figure: Detail of output map
-</center>
-
-
-<h2>EXAMPLE</h2>
-
-<div class="code"><pre>
-r.futures -s control_file=demand.txt employ_attraction=d2urbkm interchange_distance=d2intkm road_density=d2rdskm undeveloped=lc96 development_pressure=gdp cons_weight=weight_1 additional_variable_files=slope n_dev_neighbourhood=10 devpot_params=devpotParams.cfg algorithm=stochastic2 parcel_size_file=patch_sizes.txt discount_factor=0.6 probability_lookup_file=probLookup.csv sort_probs=0 patch_factor=1 patch_mean=0.4 patch_range=0.08 num_neighbors=4 seed_search=2 development_pressure_approach=2 alpha=2 scaling_factor=10 num_regions=9 index_file=index control_file_all=demand.txt output=final_results output_series=development
-</pre></div>
-
-
-<h2>KNOWN ISSUES</h2>
-
-Module holds all data in the memory, so for large areas the memory requirements
-are high.
-
-Only part of the options is actually working now. The older simpler
-implementations of the algorithms are not available.
-
-The inputs are not documented.
-
-
-<h2>REFERENCES</h2>
-
-<ul>
-<li>
- Meentemeyer, R. K., Tang, W., Dorning, M. A, Vogler, J. B., Cunniffe, N. J., Shoemaker D. A.,
- <a href="http://dx.doi.org/10.1080/00045608.2012.707591">FUTURES: Multilevel Simulations of Emerging Urban-Rural Landscape Structure Using a Stochastic Patch-Growing Algorithm</a>
- Annals of the Association of American Geographers
- Vol. 103, Iss. 4, 2013
-<li>
- Dorning, M. A, Koch J., Shoemaker D. A, Meentemeyer, R. K.,
- Simulating urbanization scenarios under various landscape planning policies reveals trade-offs between disparate conservation goals
- [in press]
-</ul>
-
-
-<h2>SEE ALSO</h2>
-
-<em>
- <a href="r.mapcalc.html">r.mapcalc</a>,
- <a href="r.reclass.html">r.reclass</a>,
- <a href="r.rescale.html">r.rescale</a>
-</em>
-
-
-<h2>AUTHORS</h2>
-
-<p>
-<em>Corresponding author:</em><br>
-Ross K. Meentemeyer, rkmeente ncsu edu,
-<a href="http://geospatial.ncsu.edu/">Center for Geospatial Analytics, NCSU</a>
-
-<p>
-<em>Original standalone version:</em><br>
-
-Ross K. Meentemeyer *<br>
-Wenwu Tang *<br>
-Monica A. Dorning *<br>
-John B. Vogler *<br>
-Nik J. Cunniffe *<br>
-Douglas A. Shoemaker *<br>
-Jennifer A. Koch **<br>
-
-<br>
-* Department of Geography and Earth Sciences, UNC Charlotte<br>
-** <a href="http://geospatial.ncsu.edu/">Center for Geospatial Analytics, NCSU</a><br>
-
-<p>
-<em>Port to GRASS GIS and GRASS-specific additions:</em><br>
-
-Vaclav Petras, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a><br>
-
-
-<p>
-<i>Last changed: $Date$</i>
Added: grass-addons/grass7/raster/r.futures/r.futures.html
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.html (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.html 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,156 @@
+<h2>DESCRIPTION</h2>
+
+<em>r.futures.*</em> is an implementation of
+FUTure Urban-Regional Environment Simulation (FUTURES)
+which is a model for multilevel simulations of emerging urban-rural
+landscape structure.
+
+The FUTure Urban-Regional Environment Simulation (FUTURES) produces
+regional projections of landscape patterns using coupled submodels
+that integrate nonstationary drivers of land change: per capita demand (DEMAND submodel),
+site suitability (POTENTIAL submodel), and the spatial structure of conversion events (PGA submodel).
+
+<h3>Submodels</h3>
+<dl>
+ <dt><em>DEMAND</em></dt>
+ <dd>DEMAND estimates the rate of per capita land consumption
+ specific to each subregion. Projections of land consumption are based
+ on extrapolations between historical changes in population
+ and land conversion based on scenarios of future population growth.
+ How to construct the per capita demand relationship for subregions depends
+ on user's preferences and data availability.
+ Land area conversion over time can be derived for the USA e.g.
+ from National Land Cover Dataset.
+ A possible implementation of the DEMAND submodel will be available as module
+ r.futures.demand.</dd>
+
+ <dt><em>POTENTIAL</em></dt>
+ <dd>The POTENTIAL submodel uses site suitability modeling approaches
+ to quantify spatial gradients of land development potential.
+ The model uses multilevel logistic regression to
+ account for hierarchical characteristics of the land use system
+ (variation among jurisdictional structures) and
+ account for divergent relationships between predictor and response variables.
+ To generate a binary, developed-undeveloped response variable
+ using a stratified-random sample,
+ see module <a href="r.sample.category.html">r.sample.category</a>.
+ The coefficients for the statistical model that are used to
+ calculate the value of development potential can be derived in different ways,
+ one possible implementation using R will be available as module r.futures.potential.
+ One of the predictor variables is development pressure (computed using
+ <a href="r.futures.devpressure.html">r.futures.devpressure</a>)
+ which is updated each step and thus creates positive feedback
+ resulting in new development attracting even more development.
+ </dd>
+
+ <dt><em>PGA</em></dt>
+ <dd>Patch-Growing Algorithm is a stochastic algorithm, which
+ simulates undeveloped to developed land change by iterative site selection
+ and a contextually aware region growing mechanism.
+ Simulations of change at each time step feed development pressure back
+ to the POTENTIAL submodel, influencing site suitability for the next step.
+ PGA is implemented in <a href="r.futures.pga.html">r.futures.pga</a>.</dd>
+</dl>
+
+<center>
+<!-- Diagram author: Monica Dorning-->
+<img width="50%" src="FUTURES_inputs_diagram.png">
+<p>
+Figure: FUTURES submodels and input data
+</center>
+
+<h3>Input data</h3>
+We need to collect the following data:
+<dl>
+ <dt><em>Study extent and resolution</em></dt>
+ <dd>Specified with <a href="g.region.html">g.region</a> command.</dd>
+ <dt><em>Subregions</em></dt>
+ <dd>FUTURES is designed to capture variation across specified subregions
+ within the full study extent. Subregions can be for example counties.
+ DEMAND and POTENTIAL can both be specified
+ according to subregions.
+ Subregion raster map contains the subregion index for each cell as integer starting from 1.
+ If you do not wish to model by subregion, all values in this map should be 1.</dd>
+
+ <dt><em>Population data</em></dt>
+ <dd>DEMAND submodel needs historical population data for each subregion
+ for reference period and population projections for the simulated period.</dd>
+ <dt><em>Development change</em></dt>
+ <dd>Based on the change in developed cells in the beginning and
+ end of the reference period, and the population data,
+ DEMAND computes how many cells to convert for each region at each time step.
+ Development change is also used for deriving the patch sizes and shape in calibration step
+ (see <a href="r.futures.calib.html">r.futures.calib</a>) to be used in PGA submodel.
+ DEMAND and PGA require a raster map representing the starting state
+ of the landscape at the beginning of the simulation (developed = 1,
+ available for development = 0, excluded from development as NULLs).</dd>
+ <dt><em>Predictors</em></dt>
+ <dd>Development potential (POTENTIAL submodel) requires
+ a set of uncorrelated predictors (raster maps) driving the land change.
+ These can include distance to roads, distance to interchanges, slope, ...</dd>
+ <dt><em>Development pressure</em></dt>
+ <dd>The development pressure variable is one of the predictors,
+ but it is recalculated at each time step to allow for positive feedback
+ (new development attracts more development). For computing development pressure,
+ see <a href="r.futures.devpressure.html">r.futures.devpressure.</a></dd>
+</dl>
+
+<p>
+<center>
+<img src="r_futures.png">
+<p>
+Figure: FUTURES simulation result
+</center>
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li>
+ Meentemeyer, R. K., Tang, W., Dorning, M. A., Vogler, J. B., Cunniffe, N. J., & Shoemaker, D. A. (2013).
+ <a href="http://dx.doi.org/10.1080/00045608.2012.707591">FUTURES: Multilevel Simulations of Emerging
+ Urban-Rural Landscape Structure Using a Stochastic Patch-Growing Algorithm</a>.
+ Annals of the Association of American Geographers, 103(4), 785-807.
+<li>Dorning, M. A., Koch, J., Shoemaker, D. A., & Meentemeyer, R. K. (2015).
+ <a href="http://dx.doi.org/10.1016/j.landurbplan.2014.11.011">Simulating urbanization scenarios reveals
+ tradeoffs between conservation planning strategies</a>.
+ Landscape and Urban Planning, 136, 28-39.</li>
+</ul>
+
+
+<h2>SEE ALSO</h2>
+<em><a href="r.futures.pga.html">r.futures.pga</a></em>,
+<em><a href="r.futures.devpressure.html">r.futures.devpressure</a></em>,
+<em><a href="r.futures.calib.html">r.futures.calib</a></em>,
+<em><a href="r.sample.category.html">r.sample.category</a></em>
+
+
+<h2>AUTHORS</h2>
+
+<p>
+<em>Corresponding author:</em><br>
+Ross K. Meentemeyer, rkmeente ncsu edu,
+<a href="http://geospatial.ncsu.edu/">Center for Geospatial Analytics, NCSU</a>
+
+<p>
+<em>Original standalone version:</em><br>
+
+Ross K. Meentemeyer *<br>
+Wenwu Tang *<br>
+Monica A. Dorning *<br>
+John B. Vogler *<br>
+Nik J. Cunniffe *<br>
+Douglas A. Shoemaker *<br>
+Jennifer A. Koch **<br>
+
+<br>
+* Department of Geography and Earth Sciences, UNC Charlotte<br>
+** <a href="http://geospatial.ncsu.edu/">Center for Geospatial Analytics, NCSU</a><br>
+
+<p>
+<em>Port to GRASS GIS and GRASS-specific additions:</em><br>
+
+Vaclav Petras, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a><br>
+
+
+<p>
+<i>Last changed: $Date: 2014-03-25 20:16:47 -0400 (Tue, 25 Mar 2014) $</i>
Copied: grass-addons/grass7/raster/r.futures/r.futures.pga/Makefile (from rev 65540, grass-addons/grass7/raster/r.futures/Makefile)
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.pga/Makefile (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.pga/Makefile 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,16 @@
+MODULE_TOPDIR = ../..
+
+PGM = r.futures.pga
+
+LIBES = $(RASTERLIB) $(GISLIB) $(MATHLIB)
+DEPENDENCIES = $(RASTERDEP) $(GISDEP)
+
+include $(MODULE_TOPDIR)/include/Make/Module.make
+
+EXTRA_CFLAGS = -Wall
+
+LINK = $(CXX)
+
+ifneq ($(strip $(CXX)),)
+default: cmd
+endif
Added: grass-addons/grass7/raster/r.futures/r.futures.pga/incentive.png
===================================================================
(Binary files differ)
Property changes on: grass-addons/grass7/raster/r.futures/r.futures.pga/incentive.png
___________________________________________________________________
Added: svn:mime-type
+ application/octet-stream
Copied: grass-addons/grass7/raster/r.futures/r.futures.pga/main.cpp (from rev 65540, grass-addons/grass7/raster/r.futures/main.cpp)
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.pga/main.cpp (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.pga/main.cpp 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,2283 @@
+
+/****************************************************************************
+ *
+ * MODULE: r.futures.pga
+ * AUTHOR(S): Ross K. Meentemeyer
+ * Wenwu Tang
+ * Monica A. Dorning
+ * John B. Vogler
+ * Nik J. Cunniffe
+ * Douglas A. Shoemaker
+ * Jennifer A. Koch
+ * Vaclav Petras <wenzeslaus gmail com>
+ * (See the manual page for details and references.)
+ *
+ * PURPOSE: Simulation of urban-rural landscape structure (FUTURES model)
+ *
+ * COPYRIGHT: (C) 2013-2014 by Meentemeyer et al.
+ *
+ * This program is free software: you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation, either version 3
+ * of the License, or (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public
+ * License along with this program. If not, see
+ * <http://www.gnu.org/licenses/> or read the file COPYING that
+ * comes with GRASS GIS for details.
+ *
+ *****************************************************************************/
+
+/**
+ \file main.c
+
+ The main file containing both the model code and the data handing part.
+
+ A lot of parts of the code are dead code. May be useful as an alternative
+ algorithms. Some of the dead code is commented out and might need to
+ be revived if other parts will show as useful.
+
+ The language of the code is subject to change. The goal is to use either
+ C or C++, not both mixed as it is now.
+
+ Major refactoring of the code is expected.
+*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <math.h>
+#include <sys/time.h>
+#include <iostream>
+#include <fstream>
+#include <vector>
+
+extern "C"
+{
+#include <grass/gis.h>
+#include <grass/raster.h>
+#include <grass/glocale.h>
+}
+
+//#include "distance.h"
+
+/** used to flag cells that are off the lattice */
+#define _CELL_OUT_OF_RANGE -2
+
+/** used to flag cells that are not in this county */
+#define _CELL_OUT_OF_COUNTY -1
+
+/** used to flag cells that are valid */
+#define _CELL_VALID 1
+
+/** need to dynamically allocate this as putting all on stack will crash most compilers */
+#define _N_MAX_DYNAMIC_BUFF_LEN (1024*1024)
+#define _N_MAX_FILENAME_LEN 1024
+
+/** Max length of input line */
+#define N_MAXREADINLEN 8192
+
+/** for working out where to start reading data */
+#define _GIS_HEADER_LENGTH 6
+
+/** strictly-speaking, should probably parse these from GIS files */
+#define _GIS_NO_DATA_STRING "-9999"
+
+/** Value used to represent NULL (NoData) in some parts of the model */
+#define _GIS_NO_DATA_INT -9999
+
+/** use this for tDeveloped if cell still undeveloped */
+#define _N_NOT_YET_DEVELOPED -1
+
+/* algorithm to use */
+
+/** deterministic model */
+#define _N_ALGORITHM_DETERMINISTIC 1
+
+/** stochastic model with downweighting of logit probabilities and of devPressure over time */
+#define _N_ALGORITHM_STOCHASTIC_I 2
+#define _N_ALGORITHM_STOCHASTIC_II 3
+
+#define _MAX_RAND_FIND_SEED_FACTOR 25
+#define maxNumAddVariables 10
+
+/** maximal number of counties allowed */
+#define MAXNUM_COUNTY 50
+#define MAX_YEARS 100
+/// maximum array size for undev cells (maximum: 1840198 for a county within 16 counties)
+#define MAX_UNDEV_SIZE 2000000
+using namespace std;
+
+
+/* Wenwu Tang */
+/// use absolute paths
+char dirName[100];
+
+/// string for temporarily storage
+char tempStr[100];
+
+
+
+typedef struct
+{
+
+ /** see #define's starting _CELL_ above */
+ int nCellType;
+
+ /** x position on the lattice; static */
+ int thisX;
+
+ /** y position on the lattice; static */
+ int thisY;
+
+ /** whether this site is still undeveloped; varies. Note if bUndeveloped = 0 the values of employAttraction etc. are not to be trusted */
+ int bUndeveloped;
+ // TODO: is the following int or double (defined as double but usually casting to int in assignment, was as int in some print)
+
+ /** development pressure; varies */
+ double devPressure;
+
+ /** stores whether this site has been considered for development, to handle "protection" */
+ int bUntouched;
+
+ /** timestep on which developed (0 for developed at start, _N_NOT_YET_DEVELOPED for not developed yet ) */
+ int tDeveloped;
+
+ /** multiplicative factor on the probabilities */
+ double consWeight;
+
+ /** additional variables */
+ double additionVariable[maxNumAddVariables];
+ int index_region;
+ float devProba;
+} t_Cell;
+
+typedef struct
+{
+
+ /** id of this cell */
+ int cellID;
+
+ /** value of the logit */
+ double logitVal;
+
+ /** whether or not cell previously considered...need to use this in the sort */
+ int bUntouched;
+
+ /** to support random pick based on their logitVal */
+ double cumulProb;
+} t_Undev;
+
+typedef struct
+{
+
+ /** array of cells (see posFromXY for how they are indexed) */
+ t_Cell *asCells;
+
+ /** number of columns in the grid */
+ int maxX;
+
+ /** number of rows in the grid */
+ int maxY;
+
+ /** total number of cells */
+ int totalCells;
+
+ /** array of information on undeveloped cells */
+ t_Undev *asUndev;
+
+ /** number of cells which have not yet converted */
+ int undevSites;
+
+ /** posterior sample from parcel size distribution */
+ int *aParcelSizes;
+
+ /** number in that sample */
+ int parcelSizes;
+ std::vector < std::vector < t_Undev > >asUndevs; //WT
+ int num_undevSites[MAXNUM_COUNTY]; //WT
+} t_Landscape;
+
+typedef struct
+{
+
+ /** size of the grids */
+ int xSize;
+
+ /** size of the grids */
+ int ySize;
+
+ /** files containing the information to read in */
+ char *developedFile;
+ char *devPressureFile;
+ char *consWeightFile;
+ char *probLookupFile;
+ int nProbLookup;
+ double *adProbLookup;
+ // parameters provided in other ways
+ // /* parameters that go into regression formula */
+ // double dIntercept; /* 0.038884 */
+ // double dEmployAttraction; /* -0.0000091946 */
+ // double dInterchangeDistance; /* 0.000042021 */
+ // double dRoadDensity; /* 0.00065813 */
+ // double dDevPressure; /* -0.026190 */
+ /* size of square used to recalculate development pressure */
+
+ int nDevNeighbourhood; /** 8 (i.e. 8 in each direction, leading to 17 x 17 = 289 in total) */
+
+ /** used in writing rasters */
+ char *dumpFile;
+ char *outputSeries;
+
+ /** 1 deterministic, 2 old stochastic, 3 new stochastic */
+ int nAlgorithm;
+
+ /** file containing parcel size information */
+ char *parcelSizeFile;
+ double discountFactor; ///< for calibrate patch size
+
+ /** these parameters only relevant for new stochastic algorithm */
+ int sortProbs;
+ double patchMean;
+ double patchRange;
+ /// 4 or 8 neighbors only
+ int numNeighbors;
+ /// 1: uniform distribution 2: based on dev. proba.
+ int seedSearch;
+ int numAddVariables;
+
+ /** parameters for additional variables */
+ double addParameters[maxNumAddVariables][MAXNUM_COUNTY];
+ char *addVariableFile[maxNumAddVariables];
+ /// 1: #occurrence; 2: gravity (need to define alpha and scaling factor); 3: kernel, need to define alpha and scaling factor
+ ///
+ /// for 2: formula-> scalingFactor/power(dist,alpha)
+ /// for 3: formula-> scalingFactor*exp(-2*dist/alpha)
+ int devPressureApproach;
+ /// scale distance-based force
+ double scalingFactor;
+ ///< constraint on distance
+ double alpha;
+
+ /** parameters that go into regression formula */
+ double dIntercepts[MAXNUM_COUNTY];
+ double dV4[MAXNUM_COUNTY];
+ int num_Regions;
+
+ /** index file to run multiple regions */
+ char *indexFile;
+
+ /// control files for all regions
+ /// development demand for multiple regions
+ char *controlFileAll;
+ int devDemand[MAX_YEARS];
+ int devDemands[MAXNUM_COUNTY][MAX_YEARS];
+ /// number of simulation steps
+ int nSteps;
+} t_Params;
+
+
+//From WT: move here to be global variables
+t_Params sParams;
+t_Landscape sLandscape;
+
+int getUnDevIndex(t_Landscape * pLandscape);
+void Shuffle4(int *mat);
+double getDistance1(double x1, double y1, double x2, double y2);
+void readDevPotParams(t_Params * pParams, char *fn);
+void readIndexData(t_Landscape * pLandscape, t_Params * pParams);
+void findAndSortProbsAll(t_Landscape * pLandscape, t_Params * pParams,
+ int step);
+void updateMap1(t_Landscape * pLandscape, t_Params * pParams, int step,
+ int regionID);
+int getUnDevIndex1(t_Landscape * pLandscape, int regionID);
+
+
+void readDevPotParams(t_Params * pParams, char *fn)
+{
+ char str[200];
+ int id;
+ double di, d4, val;
+
+ // TODO: d5 is in file but unused
+ int i, j;
+ ifstream f;
+
+ f.open(fn);
+ if (!f.is_open())
+ G_fatal_error(_("Cannot open development potential parameters file <%s>"),
+ fn);
+
+ f.getline(str, 200);
+ for (i = 0; i < pParams->num_Regions; i++) {
+ // TODO: read by lines to count the variables
+ f >> id >> di >> d4;
+ pParams->dIntercepts[i] = di;
+ pParams->dV4[i] = d4;
+ for (j = 0; j < pParams->numAddVariables; j++) {
+ f >> val;
+ pParams->addParameters[j][i] = val;
+ }
+ }
+ f.close();
+}
+
+/**
+ Generate uniform number on [0,1).
+
+ Encapsulated to allow easy replacement if necessary.
+*/
+double uniformRandom()
+{
+ int nRet;
+
+ nRet = RAND_MAX;
+ while (nRet == RAND_MAX) { /* make sure never get 1.0 */
+ nRet = rand();
+ }
+ return ((double)nRet / (double)(RAND_MAX));
+}
+
+/**
+ Seed random number generator.
+
+ Encapsulated to allow easy replacement if necessary.
+*/
+void seedRandom(struct timeval ttime)
+{
+ srand((ttime.tv_sec * 100) + (ttime.tv_usec / 100));
+ //srand((unsigned int)time(NULL));
+}
+
+
+/**
+ Work out x and y position on the grid from index into siteMapping/landscape array.
+*/
+void xyFromPos(int nPos, int *pnX, int *pnY, t_Landscape * pLandscape)
+{
+ /* integer division just gives whole number part */
+ *pnX = nPos % pLandscape->maxX;
+ *pnY = nPos / pLandscape->maxX;
+}
+
+/**
+ Work out index into asCells array from location.
+*/
+int posFromXY(int nX, int nY, t_Landscape * pLandscape)
+{
+ int nPos;
+
+ /*
+ 0 1 2 .... nX-1
+ nX nX+1 nX+2 .... 2*nX-1
+ . .
+ . .
+ . .
+ (nY-1)*nX ... nX*nY-1
+ */
+ nPos = _CELL_OUT_OF_RANGE;
+ if (nX >= 0 && nX < pLandscape->maxX && nY >= 0 && nY < pLandscape->maxY) {
+ nPos = nX * pLandscape->maxY + nY;
+ }
+ return nPos;
+}
+
+/**
+ Allocate memory to store information on cells.
+*/
+int buildLandscape(t_Landscape * pLandscape, t_Params * pParams)
+{
+ int bRet;
+
+ /* blank everything out */
+ G_verbose_message("Initialization and memory allocation...");
+ bRet = 0;
+
+ pLandscape->asCells = NULL;
+ pLandscape->asUndev = NULL;
+ pLandscape->undevSites = 0;
+ pLandscape->aParcelSizes = NULL;
+ pLandscape->parcelSizes = 0;
+
+ /* set initial values across the landscape */
+
+ pLandscape->maxX = pParams->xSize;
+ pLandscape->maxY = pParams->ySize;
+ pLandscape->totalCells = pLandscape->maxX * pLandscape->maxY;
+ int a = sizeof(t_Cell);
+ pLandscape->asCells =
+ (t_Cell *) malloc(sizeof(t_Cell) * pLandscape->totalCells);
+ if (pLandscape->asCells) {
+ if (pParams->num_Regions > 1) {
+ /* just allocate enough space for every cell to be undev */
+ pLandscape->asUndev = (t_Undev *) malloc(sizeof(t_Undev) * pLandscape->totalCells);
+ if (pLandscape->asUndev) {
+ bRet = 1;
+ }
+ }
+ else
+ bRet = 1;
+ }
+ return bRet;
+}
+
+void readData4AdditionalVariables(t_Landscape * pLandscape,
+ t_Params * pParams)
+{
+ int i;
+ int ii;
+ double val;
+
+ for (i = 0; i < pParams->numAddVariables; i++) {
+ G_verbose_message("Reading predictor variables %s...", pParams->addVariableFile[i]);
+
+ /* open raster map */
+ int fd = Rast_open_old(pParams->addVariableFile[i], "");
+ RASTER_MAP_TYPE data_type = Rast_get_map_type(fd);
+ void *buffer = Rast_allocate_buf(data_type);
+
+ ii = 0;
+ for (int row = 0; row < pParams->xSize; row++) {
+ Rast_get_row(fd, buffer, row, data_type);
+ void *ptr = buffer;
+
+ for (int col = 0; col < pParams->ySize; col++,
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(data_type))) {
+ if (data_type == DCELL_TYPE)
+ val = *(DCELL *) ptr;
+ else if (data_type == FCELL_TYPE)
+ val = *(FCELL *) ptr;
+ else {
+ val = *(CELL *) ptr;
+ }
+ if (Rast_is_null_value(ptr, data_type)) {
+ pLandscape->asCells[ii].nCellType =
+ _CELL_OUT_OF_COUNTY;
+ }
+ if (pLandscape->asCells[ii].nCellType == _CELL_VALID)
+ pLandscape->asCells[ii].additionVariable[i] = val;
+ else
+ pLandscape->asCells[ii].additionVariable[i] = 0.0;
+ ii++;
+ }
+ }
+ Rast_close(fd);
+ G_verbose_message("Done");
+ }
+}
+
+void readIndexData(t_Landscape * pLandscape, t_Params * pParams)
+{
+ int ii;
+ int val;
+
+ ii = 0;
+ int fd = Rast_open_old(pParams->indexFile, "");
+
+ /* TODO: data type should always be CELL */
+ RASTER_MAP_TYPE data_type = Rast_get_map_type(fd);
+ void *buffer = Rast_allocate_buf(data_type);
+
+ G_verbose_message("Reading subregions %s", pParams->indexFile);
+ for (int row = 0; row < pParams->xSize; row++) {
+ Rast_get_row(fd, buffer, row, data_type);
+ void *ptr = buffer;
+
+ for (int col = 0; col < pParams->ySize; col++,
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(data_type))) {
+ if (Rast_is_null_value(ptr, data_type))
+ val = _GIS_NO_DATA_INT; // assign FUTURES representation of NULL value
+ else {
+ val = *(CELL *) ptr;
+ if (val > pParams->num_Regions)
+ G_fatal_error(_("Region ID (%d) inconsistent with number of regions (%d) in map <%s>"),
+ val, pParams->num_Regions,
+ pParams->indexFile);
+ if (val < 1)
+ G_fatal_error(_("Region ID (%d) is smaller then 1 in map <%s>"),
+ val, pParams->indexFile);
+ }
+
+ pLandscape->asCells[ii].index_region = val;
+ ii++;
+ }
+ }
+ G_verbose_message("Done");
+}
+
+/**
+ Read data from GIS rasters and put in correct places in memory.
+*/
+int readData(t_Landscape * pLandscape, t_Params * pParams)
+{
+ // TODO: this function should be rewritten to be nicer (unless storing in memorey will change completely)
+ // TODO: fix null handling somewhere
+ char *szBuff;
+ char szFName[_N_MAX_FILENAME_LEN];
+ int bRet, x, y, i, j;
+ double dVal;
+
+ G_verbose_message("Reading input rasters...");
+ bRet = 0;
+ szBuff = (char *)malloc(_N_MAX_DYNAMIC_BUFF_LEN * sizeof(char));
+ if (szBuff) {
+ for (j = 0; j < 3; j++) {
+ /* workaround to skip loading constraint map so that it can be omitted in input */
+ if (j == 2) {
+ if (!pParams->consWeightFile) {
+ i = 0;
+ for (int row = 0; row < pParams->xSize; row++) {
+ for (int col = 0; col < pParams->ySize; col++) {
+ pLandscape->asCells[i].consWeight = 1;
+ i++;
+ }
+ }
+ continue;
+ }
+ }
+ switch (j) { /* get correct filename */
+ case 0:
+ strcpy(szFName, pParams->developedFile);
+ break;
+ case 1:
+ strcpy(szFName, pParams->devPressureFile);
+ break;
+ case 2:
+ strcpy(szFName, pParams->consWeightFile);
+ break;
+ default:
+ G_fatal_error("readData(): shouldn't get here");
+ break;
+ }
+ G_verbose_message("%s...", szFName);
+ /* open raster map */
+ int fd = Rast_open_old(szFName, "");
+
+ i = 0;
+
+ x = 0;
+ bRet = 1; /* can only get worse from here on in */
+ RASTER_MAP_TYPE data_type = Rast_get_map_type(fd);
+ void *buffer = Rast_allocate_buf(data_type);
+
+ for (int row = 0; row < pParams->xSize; row++) {
+ y = 0;
+ Rast_get_row(fd, buffer, row, data_type);
+ void *ptr = buffer;
+
+ for (int col = 0; col < pParams->ySize; col++,
+ ptr = G_incr_void_ptr(ptr, Rast_cell_size(data_type))) {
+ CELL iVal;
+
+ if (data_type == DCELL_TYPE)
+ dVal = *(DCELL *) ptr;
+ else if (data_type == FCELL_TYPE)
+ dVal = *(FCELL *) ptr;
+ else {
+ iVal = *(CELL *) ptr; // integer value for some rasters
+ dVal = iVal;
+ }
+ if (j == 0) {
+ /* check for NO_DATA */
+ if (Rast_is_null_value(ptr, data_type)) {
+ pLandscape->asCells[i].nCellType =
+ _CELL_OUT_OF_COUNTY;
+ }
+ else {
+ pLandscape->asCells[i].nCellType = _CELL_VALID;
+ }
+ /* and put in the correct place */
+ pLandscape->asCells[i].thisX = x;
+ pLandscape->asCells[i].thisY = y;
+ /* and set the time of development */
+ pLandscape->asCells[i].tDeveloped = _GIS_NO_DATA_INT;
+ }
+ else {
+ /* clean up missing data */
+ if (Rast_is_null_value(ptr, data_type)) {
+ dVal = 0.0;
+ }
+ }
+ // TODO: how do we know that this is set?
+ if (pLandscape->asCells[i].nCellType == _CELL_VALID) {
+ /* put data into a correct place */
+ switch (j) {
+ case 0:
+ pLandscape->asCells[i].bUndeveloped = (iVal == 0) ? 1 : 0;
+ pLandscape->asCells[i].bUntouched = 1;
+ if (pLandscape->asCells[i].bUndeveloped == 1) {
+ pLandscape->asCells[i].tDeveloped =
+ _N_NOT_YET_DEVELOPED;
+ }
+ else {
+ /* already developed at the start */
+ pLandscape->asCells[i].tDeveloped = 0;
+ }
+ break;
+ case 1:
+ pLandscape->asCells[i].devPressure = (int)dVal;
+ break;
+ case 2:
+ pLandscape->asCells[i].consWeight = dVal;
+ break;
+ default:
+ G_fatal_error("readData(): shouldn't get here");
+ break;
+ }
+ }
+ y++;
+ i++;
+ }
+
+ x++;
+ }
+ if (x == pLandscape->maxX) {
+
+ }
+ else {
+ if (bRet) {
+ G_warning("readData(): x too small");
+ bRet = 0;
+ }
+ }
+
+ if (y != pLandscape->maxY) {
+ if (bRet) {
+ G_message(_("x=%d y=%d"), x, y);
+ G_message(_(" xSize=%d ySize=%d"), pParams->xSize,
+ pParams->ySize);
+ G_message(_(" maxX=%d minX=%d"), pLandscape->maxX,
+ pLandscape->maxY);
+ G_warning("readData(): y too small");
+ bRet = 0;
+ }
+ }
+ else {
+ if (bRet && i == pLandscape->totalCells) {
+ G_verbose_message("Done");
+ }
+ else {
+
+ if (bRet) {
+ G_message(_("x=%d y=%d"), x, y);
+ G_message(_(" xSize=%d ySize=%d"), pParams->xSize,
+ pParams->ySize);
+ G_message(_(" maxX=%d minX=%d"), pLandscape->maxX,
+ pLandscape->maxY);
+ G_warning("readData(): not read in enough points");
+ bRet = 0;
+ }
+
+ }
+ }
+ Rast_close(fd);
+
+ }
+ free(szBuff);
+ }
+
+ return bRet;
+}
+
+
+/**
+ Compare two t_Undev structures.
+
+ To be used in qsort().
+
+ \note This sort function is in ascending order
+ (i.e. opposite to what we normally want).
+
+ \see undevCmpReverse()
+*/
+static int undevCmp(const void *a, const void *b)
+{
+ double da = ((t_Undev *) a)->logitVal;
+ double db = ((t_Undev *) b)->logitVal;
+
+ /* make sure untouched cells preferentially go to top of sort when there is asymmetry between a and b */
+ if (((t_Undev *) a)->bUntouched && !((t_Undev *) b)->bUntouched)
+ return -1;
+ if (((t_Undev *) b)->bUntouched && !((t_Undev *) a)->bUntouched)
+ return 1;
+ /* then sort by value */
+ if (da > db)
+ return 1;
+ if (da < db)
+ return -1;
+ return 0;
+}
+
+/**
+ Compare two t_Undev structures.
+
+ This is the correct way around sorting with undevCmp() and qsort().
+ */
+static int undevCmpReverse(const void *a, const void *b)
+{
+ return -1 * undevCmp(a, b);
+}
+
+
+/**
+ Write current state of developed areas.
+
+ Called at end and dumps tDeveloped for all valid cells (NULL for all others)
+
+ \param undevelopedAsNull Represent undeveloped areas as NULLs instead of -1
+ \param developmentAsOne Represent all developed areas as 1 instead of number
+ representing the step when are was developed
+ */
+// TODO: add timestamp to maps
+void outputDevRasterStep(t_Landscape * pLandscape, t_Params * pParams,
+ const char *rasterNameBase, bool undevelopedAsNull,
+ bool developmentAsOne)
+{
+ int out_fd = Rast_open_new(rasterNameBase, CELL_TYPE);
+ CELL *out_row = Rast_allocate_c_buf();
+
+ for (int x = 0; x < pLandscape->maxX; x++) {
+ Rast_set_c_null_value(out_row, pLandscape->maxY);
+ for (int y = 0; y < pLandscape->maxY; y++) {
+ // we don't check return values since we use landscape directly
+ int cellId = posFromXY(x, y, pLandscape);
+
+ // this handles NULLs on input
+ if (pLandscape->asCells[cellId].nCellType == _CELL_VALID) {
+ CELL value = pLandscape->asCells[cellId].tDeveloped;
+
+ // this handles undeveloped cells
+ if (undevelopedAsNull && value == -1)
+ continue;
+ // this handles developed cells
+ if (developmentAsOne)
+ value = 1;
+ out_row[y] = value;
+ }
+ }
+ Rast_put_c_row(out_fd, out_row);
+ }
+ Rast_close(out_fd);
+
+ struct Colors colors;
+
+ Rast_init_colors(&colors);
+ CELL val1 = 0;
+ // TODO: the map max is 36 for 36 steps, it is correct?
+ CELL val2 = pParams->nSteps;
+
+ if (developmentAsOne) {
+ val1 = 1;
+ val2 = 1;
+ Rast_add_c_color_rule(&val1, 255, 100, 50, &val2, 255, 100, 50,
+ &colors);
+ }
+ else {
+ Rast_add_c_color_rule(&val1, 255, 100, 50, &val2, 255, 255, 0,
+ &colors);
+ }
+ if (!undevelopedAsNull) {
+ val1 = -1;
+ val2 = -1;
+ Rast_add_c_color_rule(&val1, 100, 255, 100, &val2, 100, 255, 100,
+ &colors);
+ }
+
+ const char *mapset = G_find_file2("cell", rasterNameBase, "");
+
+ if (mapset == NULL)
+ G_fatal_error(_("Raster map <%s> not found"), rasterNameBase);
+
+ Rast_write_colors(rasterNameBase, mapset, &colors);
+ Rast_free_colors(&colors);
+
+ struct History hist;
+
+ Rast_short_history(rasterNameBase, "raster", &hist);
+ Rast_command_history(&hist);
+ // TODO: store also random seed value (need to get it here, global? in Params?)
+ Rast_write_history(rasterNameBase, &hist);
+
+ G_message(_("Raster map <%s> created"), rasterNameBase);
+}
+
+char *mapNameForStep(const char *basename, const int step, const int maxSteps)
+{
+ int digits = log10(maxSteps) + 1;
+
+ return G_generate_basename(basename, step, digits, 0);
+}
+
+/**
+ Work out how many cells to transition based on the control file.
+*/
+int parseControlLine(char *szBuff, char *szStepLabel, int *pnToConvert)
+{
+ char *pPtr;
+ int i;
+ int bRet;
+
+ bRet = 0;
+ /* strip newline */
+ pPtr = strpbrk(szBuff, "\r\n");
+ if (pPtr) {
+ pPtr[0] = '\0';
+ }
+ *pnToConvert = 0;
+ i = 0;
+ pPtr = strtok(szBuff, " \t");
+ while (pPtr && i < 2) {
+ switch (i) {
+ case 0:
+ strcpy(szStepLabel, pPtr);
+ break;
+ case 1:
+ *pnToConvert = atoi(pPtr);
+ break;
+ }
+ pPtr = strtok(NULL, " \t");
+ i++;
+ }
+ /* only process this line of the control file if it is well formed */
+ if (i == 2) {
+ bRet = 1;
+ }
+ return bRet;
+}
+
+/**
+ Calculate probability of development.
+*/
+
+double getDevProbability(t_Cell * pThis, t_Params * pParams)
+{
+ double probAdd;
+ int i;
+ int id = pThis->index_region;
+
+ if (id == -9999)
+ return 0;
+ id = id - 1;
+ probAdd = pParams->dIntercepts[id];
+ //cout<<"intercept\t"<<probAdd<<endl;
+ probAdd += pParams->dV4[id] * pThis->devPressure;
+ //cout<<"devPressure: "<<pParams->dV4[id]<<endl;
+ //cout<<"devPressure: "<<pThis->devPressure<<endl;
+ //cout<<"devPressure: "<<pParams->dV4[id] * pThis->devPressure<<endl;
+ for (i = 0; i < pParams->numAddVariables; i++) {
+ probAdd += pParams->addParameters[i][id] * pThis->additionVariable[i];
+ //cout<<"additionVariable: "<<i<<"\t"<<pParams->addParameters[i][id]<<endl;
+ //cout<<"additionVariable: "<<i<<"\t"<<pThis->additionVariable[i]<<endl;
+ //cout<<"additionVariable: "<<i<<"\t"<<pParams->addParameters[i][id] * pThis->additionVariable[i]<<endl;
+ }
+ probAdd = 1.0 / (1.0 + exp(-probAdd));
+ //cout<<"The probability:";
+ //cout<<probAdd<<endl;
+ return probAdd;
+}
+
+/*
+ double getDevProbability(t_Cell *pThis, t_Params *pParams){
+ float probAdd; int i;
+ int id=pThis->index_region;
+ if(id==-9999)
+ return 0;
+ id=id-1;
+ probAdd = pParams->dIntercepts[id];
+ probAdd += pParams->dV1[id] * pThis->employAttraction;
+ probAdd += pParams->dV2[id] * pThis->interchangeDistance;
+ probAdd += pParams->dV3[id] * pThis->roadDensity;
+ probAdd += pParams->dV4[id] * pThis->devPressure;
+ for(i=0;i<2;i++){
+ probAdd+=pParams->addParameters[i][id]*pThis->additionVariable[i];
+ }
+ if (0<pThis->additionVariable[2]<=400)
+ {
+ probAdd = probAdd * 0.0005 ;
+ }
+ else {
+ if(400<pThis->additionVariable[2]<=600)
+ {probAdd = probAdd * 0.024;}
+ else {
+ if(600<pThis->additionVariable[2]<=800)
+ {probAdd = probAdd * 0.912 ;}
+ else{
+ if(800<pThis->additionVariable[2]<=1000)
+ {probAdd = probAdd * 0.055 ;}
+ else{
+ if(1000<pThis->additionVariable[2]<=1200)
+ {probAdd = probAdd * 0.006 ;}
+ else{
+ if(1200<pThis->additionVariable[2]<=1400)
+ {probAdd = probAdd * 0.002 ;}
+ else{
+
+ {probAdd = probAdd * 0.0005;}
+
+ }
+ }
+ }
+
+ }
+
+ }
+
+
+ }
+
+
+ probAdd = 1.0/(1.0 + exp(-probAdd));
+ return probAdd;
+ }
+
+ */
+
+
+/**
+ Recalculate probabilities for each unconverted cell.
+
+ Called each timestep.
+*/
+void findAndSortProbs(t_Landscape * pLandscape, t_Params * pParams,
+ int nToConvert)
+{
+ int i, lookupPos;
+ t_Cell *pThis;
+
+ /* update calcs */
+ G_verbose_message("Recalculating probabilities...");
+ pLandscape->undevSites = 0;
+ for (i = 0; i < pLandscape->totalCells; i++) {
+ pThis = &(pLandscape->asCells[i]);
+ if (pThis->nCellType == _CELL_VALID) {
+ if (pThis->bUndeveloped) {
+ if (pThis->consWeight > 0.0) {
+ /* note that are no longer just storing the logit value, but instead the probability (allows consWeight to affect sort order) */
+ pLandscape->asUndev[pLandscape->undevSites].cellID = i;
+ pLandscape->asUndev[pLandscape->undevSites].logitVal =
+ getDevProbability(pThis, pParams);
+ /* lookup table of probabilities is applied before consWeight */
+ if (pParams->nAlgorithm == _N_ALGORITHM_STOCHASTIC_II) {
+ /* replace with value from lookup table */
+ lookupPos =
+ (int)(pLandscape->asUndev[pLandscape->undevSites].
+ logitVal * (pParams->nProbLookup - 1));
+ pLandscape->asUndev[pLandscape->undevSites].logitVal =
+ pParams->adProbLookup[lookupPos];
+ // fprintf(stdout, "%f %d %f\n", pLandscape->asUndev[pLandscape->undevSites].logitVal, lookupPos, pParams->adProbLookup[lookupPos]);
+ }
+ // multiply consweight
+ pLandscape->asUndev[pLandscape->undevSites].logitVal *=
+ pThis->consWeight;
+ pLandscape->asUndev[pLandscape->undevSites].bUntouched = pThis->bUntouched; /* need to store this to put correct elements near top of list */
+ if (pLandscape->asUndev[pLandscape->undevSites].logitVal >
+ 0.0) {
+ /* only add one more to the list if have a positive probability */
+ pLandscape->undevSites++;
+ }
+ }
+ }
+ }
+ }
+ /* sort */
+ /* can only be zero for algorithm=stochastic_ii */
+ if (pParams->sortProbs) {
+ G_verbose_message("Sorting %d unconserved undeveloped sites",
+ pLandscape->undevSites);
+ qsort(pLandscape->asUndev, pLandscape->undevSites, sizeof(t_Undev),
+ undevCmpReverse);
+ }
+ else {
+ G_verbose_message("Skipping sort as choosing cells randomly");
+ }
+ //calculate cumulative probability
+ // From Wenwu Tang
+ double sum = pLandscape->asUndev[0].logitVal;
+
+ for (i = 1; i < pLandscape->undevSites; i++) {
+ pLandscape->asUndev[i].cumulProb =
+ pLandscape->asUndev[i - 1].cumulProb +
+ pLandscape->asUndev[i].logitVal;
+ sum = sum + pLandscape->asUndev[i].logitVal;
+ }
+ for (i = 0; i < pLandscape->undevSites; i++) {
+ pLandscape->asUndev[i].cumulProb =
+ pLandscape->asUndev[i].cumulProb / sum;
+ }
+}
+
+
+/**
+ Helper structur for neighbour grabbing algorithm
+ */
+typedef struct
+{
+ double probAdd;
+ int cellID;
+ int newInList;
+ // Wenwu Tang
+ /// distance to the center
+ double distance;
+} t_candidateNeighbour;
+
+/**
+ Helper structur for neighbour grabbing algorithm
+*/
+typedef struct
+{
+ double maxProb;
+ int nCandidates;
+ int nSpace;
+ t_candidateNeighbour *aCandidates;
+} t_neighbourList;
+
+#define _N_NEIGHBOUR_LIST_BLOCK_SIZE 20
+
+int addNeighbourIfPoss(int x, int y, t_Landscape * pLandscape,
+ t_neighbourList * pNeighbours, t_Params * pParams)
+{
+ int i, thisPos, mustAdd, listChanged, lookupPos;
+ double probAdd;
+ t_Cell *pThis;
+
+ listChanged = 0;
+ thisPos = posFromXY(x, y, pLandscape);
+ if (thisPos != _CELL_OUT_OF_RANGE) {
+ pThis = &(pLandscape->asCells[thisPos]);
+ if (pThis->nCellType == _CELL_VALID) {
+ pThis->bUntouched = 0;
+ if (pThis->bUndeveloped) {
+ if (pThis->consWeight > 0.0) {
+ /* need to add this cell... */
+
+ /* ...either refresh its element in list if already there */
+ mustAdd = 1;
+ for (i = 0; mustAdd && i < pNeighbours->nCandidates; i++) {
+ if (pNeighbours->aCandidates[i].cellID == thisPos) {
+ pNeighbours->aCandidates[i].newInList = 1;
+ mustAdd = 0;
+ listChanged = 1;
+ }
+ }
+ /* or add it on the end, allocating space if necessary */
+ if (mustAdd) {
+ if (pNeighbours->nCandidates == pNeighbours->nSpace) {
+ pNeighbours->nSpace +=
+ _N_NEIGHBOUR_LIST_BLOCK_SIZE;
+ pNeighbours->aCandidates =
+ (t_candidateNeighbour *)
+ realloc(pNeighbours->aCandidates,
+ pNeighbours->nSpace *
+ sizeof(t_candidateNeighbour));
+ if (!pNeighbours->aCandidates) {
+ G_fatal_error("Memory error in addNeighbourIfPoss()");
+ }
+ }
+ pNeighbours->aCandidates[pNeighbours->nCandidates].
+ cellID = thisPos;
+ pNeighbours->aCandidates[pNeighbours->nCandidates].
+ newInList = 1;
+
+ /* note duplication of effort in here recalculating the probabilities, but didn't store them in an accessible way */
+ /*
+ probAdd = pParams->dIntercept;
+ probAdd += pParams->dDevPressure * pThis->devPressure;
+ probAdd += pParams->dEmployAttraction * pThis->employAttraction;
+ probAdd += pParams->dInterchangeDistance * pThis->interchangeDistance;
+ probAdd += pParams->dRoadDensity * pThis->roadDensity;
+ probAdd = 1.0/(1.0 + exp(probAdd));
+ */
+ probAdd = getDevProbability(pThis, pParams);
+ /* replace with value from lookup table */
+ lookupPos =
+ (int)(probAdd * (pParams->nProbLookup - 1));
+ probAdd = pParams->adProbLookup[lookupPos];
+ probAdd *= pThis->consWeight;
+ pNeighbours->aCandidates[pNeighbours->nCandidates].
+ probAdd = probAdd;
+ /* only actually add it if will ever transition */
+ if (probAdd > 0.0) {
+ pNeighbours->nCandidates++;
+ if (probAdd > pNeighbours->maxProb) {
+ pNeighbours->maxProb = probAdd;
+ }
+ listChanged = 1;
+ }
+ }
+ }
+ }
+ }
+ }
+ return listChanged;
+}
+
+// From Wenwu
+
+/**
+ Sort according to the new flag and conversion probability.
+*/
+static int sortNeighbours(const void *pOne, const void *pTwo)
+{
+ t_candidateNeighbour *pNOne = (t_candidateNeighbour *) pOne;
+ t_candidateNeighbour *pNTwo = (t_candidateNeighbour *) pTwo;
+
+ /*
+ if(pNOne->newInList > pNTwo->newInList)
+ {
+ return -1;
+ }
+ if(pNTwo->newInList > pNOne->newInList)
+ {
+ return 1;
+ }
+ */
+ float p = rand() / (float)RAND_MAX;
+
+ /*
+ if(p<0.0){
+ if(pNOne->distance> pNTwo->distance)
+ {
+ return 1;
+ }
+ if(pNTwo->distance > pNOne->distance)
+ {
+ return -1;
+ }
+ }
+ */
+ float alpha = 0.5; // assign this from parameter file
+
+ alpha = (sParams.patchMean) - (sParams.patchRange) * 0.5;
+ alpha += p * sParams.patchRange;
+ //alpha=0;
+ float p1, p2;
+
+ p1 = pNOne->probAdd / pow(pNOne->distance, alpha);
+ p2 = pNTwo->probAdd / pow(pNTwo->distance, alpha);
+
+ //p1=alpha*pNOne->probAdd+(1-alpha)*pow(pNOne->distance,-2);
+ //p2=alpha*pNTwo->probAdd+(1-alpha)*pow(pNTwo->distance,-2);
+
+ if (p1 > p2) {
+ return -1;
+ }
+ if (p2 > p1) {
+ return 1;
+ }
+ return 0;
+}
+
+/**
+ randomly shuffle the 4 neighbors
+*/
+void Shuffle4(int *mat)
+{
+ int proba[4], flag[4];
+ int size = 4, i;
+
+ for (i = 0; i < size; i++) {
+ proba[i] = rand() % 1000;
+ flag[i] = 1;
+ }
+
+ int numChecked = 0;
+ int max;
+ int index = 0;
+
+ while (numChecked != 4) {
+ max = -9999;
+ for (i = 0; i < size; i++) {
+ if (flag[i] != 0) {
+ if (proba[i] > max) {
+ max = proba[i];
+ index = i;
+ }
+ }
+ }
+ flag[index] = 0;
+ mat[numChecked] = index;
+ numChecked++;
+ }
+}
+
+void findNeighbours(int nThisID, t_Landscape * pLandscape,
+ t_neighbourList * pNeighbours, t_Params * pParams)
+{
+ t_Cell *pThis;
+ int listChanged = 0;
+
+ // int idmat[4]; idmat[0]=0;idmat[1]=1;idmat[2]=2;idmat[3]=3;
+ /* try to add the neighbours of this one */
+ pThis = &(pLandscape->asCells[nThisID]);
+ //note: if sorted, then shuffle is no need
+ /*
+ Shuffle4(&idmat[0]);
+ for(i=0;i<4;i++){
+ if(idmat[i]==0)
+ listChanged+=addNeighbourIfPoss(pThis->thisX-1,pThis->thisY,pLandscape,pNeighbours,pParams);
+ else if(idmat[i]==1)
+ listChanged+=addNeighbourIfPoss(pThis->thisX+1,pThis->thisY,pLandscape,pNeighbours,pParams);
+ else if(idmat[i]==2)
+ listChanged+=addNeighbourIfPoss(pThis->thisX,pThis->thisY-1,pLandscape,pNeighbours,pParams);
+ else if (idmat[i]==3)
+ listChanged+=addNeighbourIfPoss(pThis->thisX,pThis->thisY+1,pLandscape,pNeighbours,pParams);
+ }
+ */
+
+ listChanged += addNeighbourIfPoss(pThis->thisX - 1, pThis->thisY, pLandscape, pNeighbours, pParams); // left
+ listChanged += addNeighbourIfPoss(pThis->thisX + 1, pThis->thisY, pLandscape, pNeighbours, pParams); // right
+ listChanged += addNeighbourIfPoss(pThis->thisX, pThis->thisY - 1, pLandscape, pNeighbours, pParams); // down
+ listChanged += addNeighbourIfPoss(pThis->thisX, pThis->thisY + 1, pLandscape, pNeighbours, pParams); // up
+ if (sParams.numNeighbors == 8) {
+ listChanged +=
+ addNeighbourIfPoss(pThis->thisX - 1, pThis->thisY - 1, pLandscape,
+ pNeighbours, pParams);
+ listChanged +=
+ addNeighbourIfPoss(pThis->thisX - 1, pThis->thisY + 1, pLandscape,
+ pNeighbours, pParams);
+ listChanged +=
+ addNeighbourIfPoss(pThis->thisX + 1, pThis->thisY - 1, pLandscape,
+ pNeighbours, pParams);
+ listChanged +=
+ addNeighbourIfPoss(pThis->thisX + 1, pThis->thisY + 1, pLandscape,
+ pNeighbours, pParams);
+ }
+ /* if anything has been altered then resort */
+
+ if (listChanged) {
+ //qsort(pNeighbours->aCandidates,pNeighbours->nCandidates,sizeof(t_candidateNeighbour),sortNeighbours);
+ }
+}
+
+double getDistance(int id1, int id2, t_Landscape * pLandscape)
+{
+ double x1, x2, y1, y2, result = 0;
+
+ x1 = pLandscape->asCells[id1].thisX;
+ x2 = pLandscape->asCells[id2].thisX;
+ y1 = pLandscape->asCells[id1].thisY;
+ y2 = pLandscape->asCells[id2].thisY;
+ result = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
+ return result;
+}
+
+double getDistance1(double x1, double y1, double x2, double y2)
+{
+ double result = 0;
+
+ result = sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2));
+ return result;
+}
+
+int newPatchFinder(int nThisID, t_Landscape * pLandscape, t_Params * pParams,
+ int *anToConvert, int nWantToConvert)
+{
+ int bTrying, i, nFound, j;
+ t_neighbourList sNeighbours;
+
+ memset(&sNeighbours, 0, sizeof(t_neighbourList));
+ anToConvert[0] = nThisID;
+ pLandscape->asCells[nThisID].bUntouched = 0;
+ pLandscape->asCells[nThisID].bUndeveloped = 0;
+ nFound = 1;
+ findNeighbours(nThisID, pLandscape, &sNeighbours, pParams);
+ for (i = 0; i < sNeighbours.nCandidates; i++) {
+ sNeighbours.aCandidates[i].distance =
+ getDistance(nThisID, sNeighbours.aCandidates[i].cellID,
+ pLandscape);
+ }
+
+ while (nFound < nWantToConvert && sNeighbours.nCandidates > 0) {
+ i = 0;
+ bTrying = 1;
+ while (bTrying) {
+ if (uniformRandom() < sNeighbours.aCandidates[i].probAdd) {
+ anToConvert[nFound] = sNeighbours.aCandidates[i].cellID;
+ /* flag that it is about to develop */
+ pLandscape->asCells[anToConvert[nFound]].bUndeveloped = 0;
+ /* remove this one from the list by copying down everything above it */
+ for (j = i + 1; j < sNeighbours.nCandidates; j++) {
+ sNeighbours.aCandidates[j - 1].cellID =
+ sNeighbours.aCandidates[j].cellID;
+ sNeighbours.aCandidates[j - 1].newInList =
+ sNeighbours.aCandidates[j].newInList;
+ sNeighbours.aCandidates[j - 1].probAdd =
+ sNeighbours.aCandidates[j].probAdd;
+ }
+ /* reduce the size of the list */
+ sNeighbours.nCandidates--;
+ sNeighbours.maxProb = 0.0;
+ for (j = 0; j < sNeighbours.nCandidates; j++) {
+ if (sNeighbours.aCandidates[j].probAdd >
+ sNeighbours.maxProb) {
+ sNeighbours.maxProb =
+ sNeighbours.aCandidates[j].probAdd;
+ }
+ }
+ /* add its neighbours */
+ findNeighbours(anToConvert[nFound], pLandscape, &sNeighbours,
+ pParams);
+ nFound++;
+ bTrying = 0;
+ }
+ else {
+ sNeighbours.aCandidates[i].newInList = 0;
+ }
+ if (bTrying) {
+ i++;
+ if (i == sNeighbours.nCandidates) {
+ i = 0;
+ }
+ }
+ else {
+ /* always resort the list if have just added, to let new elements bubble to the top */
+ if (sNeighbours.nCandidates > 0) {
+ // int z;
+ for (i = 0; i < sNeighbours.nCandidates; i++) {
+ sNeighbours.aCandidates[i].distance =
+ getDistance(nThisID,
+ sNeighbours.aCandidates[i].cellID,
+ pLandscape);
+ }
+
+ qsort(sNeighbours.aCandidates, sNeighbours.nCandidates,
+ sizeof(t_candidateNeighbour), sortNeighbours);
+#if 0
+ for (z = 0; z < sNeighbours.nCandidates; z++) {
+ fprintf(stdout, "%d %d %f\n", z,
+ sNeighbours.aCandidates[z].newInList,
+ sNeighbours.aCandidates[z].probAdd);
+ }
+ fprintf(stdout, "***\n");
+#endif
+ }
+ }
+ }
+ }
+ if (sNeighbours.nSpace) { /* free any allocated memory */
+ free(sNeighbours.aCandidates);
+ }
+ // fprintf(stdout, "looking for %d found %d\n", nWantToConvert,nFound);
+ return nFound;
+}
+
+int convertCells(t_Landscape * pLandscape, t_Params * pParams, int nThisID,
+ int nStep, int bAllowTouched, int bDeterministic)
+{
+ int i, x, y, nCell, nDone, nToConvert, nWantToConvert;
+ int *anToConvert;
+ t_Cell *pThis, *pThat;
+
+ nDone = 0;
+ /* in here need to build list of near neigbours to loop over */
+ nWantToConvert = pLandscape->aParcelSizes[(int)
+ (uniformRandom() *
+ pLandscape->parcelSizes)];
+ anToConvert = (int *)malloc(sizeof(int) * nWantToConvert);
+ if (anToConvert) {
+ /* in here goes code to fill up list of neighbours */
+ nToConvert =
+ newPatchFinder(nThisID, pLandscape, pParams, anToConvert,
+ nWantToConvert);
+ /* will actually always be the case */
+ if (nToConvert > 0) {
+ // fprintf(stdout, "wanted %d, found %d\n", nWantToConvert, nToConvert);
+ for (i = 0; i < nToConvert; i++) {
+ /* convert, updating dev press on neighbours */
+ pThis = &(pLandscape->asCells[anToConvert[i]]);
+ pThis->bUndeveloped = 0;
+ pThis->bUntouched = 0;
+ pThis->tDeveloped = nStep + 1;
+ nDone++;
+ // fprintf(stdout, "%d total=%f dev=%f employment=%f interchange=%f roads=%f\n", pLandscape->asUndev[i].cellID, pLandscape->asUndev[i].logitVal, pParams->dDevPressure * pThis->devPressure, pParams->dEmployAttraction * pThis->employAttraction, pParams->dInterchangeDistance * pThis->interchangeDistance, pParams->dRoadDensity * pThis->roadDensity);
+ float dist = 0;
+ float force = 0;
+
+ for (x = pThis->thisX - pParams->nDevNeighbourhood;
+ x <= pThis->thisX + pParams->nDevNeighbourhood; x++) {
+ for (y = pThis->thisY - pParams->nDevNeighbourhood;
+ y <= pThis->thisY + pParams->nDevNeighbourhood;
+ y++) {
+ nCell = posFromXY(x, y, pLandscape);
+ if (nCell != _CELL_OUT_OF_RANGE) {
+ pThat = &(pLandscape->asCells[nCell]);
+ dist =
+ getDistance1(pThis->thisX, pThis->thisY, x,
+ y);
+ if (dist > pParams->nDevNeighbourhood)
+ continue;
+ force = 0;
+ if (pParams->devPressureApproach == 1)
+ force = 1; // development pressure increases by 1
+ else if (pParams->devPressureApproach == 2) {
+ if (dist > 0)
+ force =
+ pParams->scalingFactor / pow(dist,
+ pParams->
+ alpha);
+ }
+ else {
+ force =
+ pParams->scalingFactor * exp(-2 * dist /
+ pParams->
+ alpha);
+ }
+ pThat->devPressure = pThat->devPressure + force;
+ }
+ }
+ }
+ }
+ }
+ free(anToConvert);
+ }
+ return nDone;
+}
+
+
+void readDevDemand(t_Params * pParams)
+{
+ ifstream f1;
+
+ /*
+ f1.open(pParams->controlFile);
+ int counter=0;
+ int val1,val2;
+ while(!f1.eof()){
+ f1>>val1>>val2;
+ pParams->devDemand[counter]=val2;
+ counter++;
+ }
+ pParams->nSteps=counter;
+ f1.close();
+ */
+ //read land demand for each region
+ int i;
+ int val1;
+ char str[200];
+
+ f1.open(pParams->controlFileAll);
+ f1 >> str >> val1;
+ pParams->nSteps = val1;
+ f1.getline(str, 200);
+ f1.getline(str, 200); //skip the header
+ int ii;
+
+ for (ii = 0; ii < pParams->nSteps; ii++) {
+ f1 >> val1;
+ for (i = 0; i < pParams->num_Regions; i++) {
+ f1 >> val1;
+ pParams->devDemands[i][ii] = val1;
+ }
+
+ }
+ f1.close();
+
+}
+
+void initializeUnDev(t_Landscape * pLandscape, t_Params * pParams)
+{
+ int i;
+
+ pLandscape->asUndevs.clear();
+ pLandscape->asUndevs.resize(pParams->num_Regions,
+ std::vector < t_Undev > (MAX_UNDEV_SIZE));
+ for (i = 0; i < pParams->num_Regions; i++) {
+ pLandscape->num_undevSites[i] = 0;
+ }
+}
+
+void finalizeUnDev(t_Landscape * pLandscape, t_Params * pParams)
+{
+
+}
+
+/*
+ main routine to run models on multiple regions // WTang
+ */
+void updateMapAll(t_Landscape * pLandscape, t_Params * pParams)
+{
+ //readDevDemand(pParams);
+ int i, j;
+
+ //initialize undev arrays
+ initializeUnDev(pLandscape, pParams);
+
+ //iterate each step (year)
+ for (i = 0; i < pParams->nSteps; i++) {
+ G_verbose_message("Processing step %d from %d", i + 1, pParams->nSteps);
+ // for each sub-region, find and update conversion probability (conservation weight applied)
+ findAndSortProbsAll(pLandscape, pParams, i);
+ for (j = 0; j < pParams->num_Regions; j++)
+ //j=1;
+ updateMap1(pLandscape, pParams, i, j);
+ if (pParams->outputSeries)
+ outputDevRasterStep(pLandscape, pParams,
+ mapNameForStep(pParams->outputSeries, i,
+ pParams->nSteps), true, true);
+ }
+ outputDevRasterStep(pLandscape, pParams, pParams->dumpFile, false, false);
+ finalizeUnDev(pLandscape, pParams);
+}
+
+void updateMap1(t_Landscape * pLandscape, t_Params * pParams, int step,
+ int regionID)
+{
+ int i, nToConvert, nStep, nDone, bAllowTouched, nRandTries, nExtra;
+ double dProb;
+ t_Cell *pThis;
+
+ nExtra = 0;
+
+ nStep = step;
+
+ //fprintf(stdout, "entered updateMap()\n");
+ nToConvert = pParams->devDemands[regionID][step];
+ //fprintf(stdout, "\tdoing step %s...controlFile requests conversion of %d cells\n", szStepLabel, nToConvert);
+ if (nExtra > 0) {
+ if (nToConvert - nExtra > 0) {
+ nToConvert -= nExtra;
+ nExtra = 0;
+ }
+ else {
+ nExtra -= nToConvert;
+ nToConvert = 0;
+ }
+ }
+ G_debug(1, "After accounting for extra cells, attempt %d cells", nToConvert);
+ /* if have cells to convert this step */
+ if (nToConvert > 0) {
+ //findAndSortProbs(pLandscape, pParams, nToConvert);
+ /* if not enough cells to convert then alter number required */
+ if (nToConvert > pLandscape->num_undevSites[regionID]) {
+ G_warning("Not enough undeveloped sites... converting all");
+ nToConvert = pLandscape->num_undevSites[regionID];
+ }
+ /* update either in deterministic or stochastic fashion */
+ //fprintf(stdout, "\t\tupdating map\n");
+ switch (pParams->nAlgorithm) {
+ case _N_ALGORITHM_DETERMINISTIC: /* deterministic */
+ nDone = 0;
+ for (i = 0; i < nToConvert && nDone < nToConvert; i++) {
+ nDone +=
+ convertCells(pLandscape, pParams,
+ pLandscape->asUndev[i].cellID, nStep, 1, 1);
+ }
+ break;
+ case _N_ALGORITHM_STOCHASTIC_II: /* stochastic */
+ nDone = 0;
+ i = 0;
+ nRandTries = 0;
+ bAllowTouched = 0;
+ /* loop until done enough cells...might need multiple passes */
+ while (nDone < nToConvert) {
+ if (i == pLandscape->num_undevSites[regionID]) {
+ /* if at the end of the grid, just loop around again until done */
+ i = 0;
+ /* allow previously considered cells if you have to */
+ bAllowTouched = 1;
+ }
+ /* if tried too many randomly in this step, give up on idea of only letting untouched cells convert */
+ if (nRandTries > _MAX_RAND_FIND_SEED_FACTOR * nToConvert) {
+ bAllowTouched = 1;
+ }
+ if (pParams->sortProbs) {
+ /* if sorted then choose the top cell and do nothing */
+ }
+ else {
+ /* otherwise give a random undeveloped cell a go */
+ if (sParams.seedSearch == 1)
+ i = (int)(uniformRandom() *
+ pLandscape->num_undevSites[regionID]);
+ //pick one according to their probability
+ else
+ G_debug(3, "nDone=%d, toConvert=%d", nDone,
+ nToConvert);
+ i = getUnDevIndex1(pLandscape, regionID);
+ }
+ pThis =
+ &(pLandscape->asCells
+ [pLandscape->asUndevs[regionID][i].cellID]);
+
+ /* need to check is still undeveloped */
+ if (pThis->bUndeveloped) {
+ if (bAllowTouched || pThis->bUntouched) {
+ /* Doug's "back to front" logit */
+ // dProb = 1.0/(1.0 + exp(pLandscape->asUndev[i].logitVal));
+ dProb = pLandscape->asUndevs[regionID][i].logitVal;
+ if (uniformRandom() < dProb) {
+ nDone +=
+ convertCells(pLandscape, pParams,
+ pLandscape->asUndevs[regionID]
+ [i].cellID, nStep, bAllowTouched,
+ 0);
+ }
+ }
+ pThis->bUntouched = 0;
+ }
+ if (pParams->sortProbs) {
+ i++;
+ }
+ else {
+ nRandTries++;
+ }
+ }
+ break;
+ default:
+ G_fatal_error("Unknown algorithm...exiting");
+ break;
+ }
+ G_debug(1, "Converted %d sites", nDone);
+ nExtra += (nDone - nToConvert);
+ G_debug(1, "%d extra sites knocked off next timestep", nExtra);
+ }
+}
+
+
+/**
+ Check Doug's calculation of devPressure.
+
+ No longer called.
+ */
+void testDevPressure(t_Landscape * pLandscape, t_Params * pParams)
+{
+ int x, y, xN, yN, cellID, nCell, testPressure;
+
+ for (y = 0; y < pLandscape->maxY; y++) {
+ for (x = 0; x < pLandscape->maxX; x++) {
+ cellID = posFromXY(x, y, pLandscape);
+ if (cellID != _CELL_OUT_OF_RANGE) { /* should never happen */
+ if (pLandscape->asCells[cellID].nCellType == _CELL_VALID) {
+ if (pLandscape->asCells[cellID].bUndeveloped) {
+ if (pLandscape->asCells[cellID].devPressure) {
+ fprintf(stdout,
+ "(x,y)=(%d,%d) cellType=%d bUndev=%d devPress=%f\n",
+ x, y,
+ pLandscape->asCells[cellID].nCellType,
+ pLandscape->asCells[cellID].bUndeveloped,
+ pLandscape->asCells[cellID].devPressure);
+ testPressure = 0;
+ for (xN = x - pParams->nDevNeighbourhood;
+ xN <= x + pParams->nDevNeighbourhood; xN++) {
+ for (yN = y - pParams->nDevNeighbourhood;
+ yN <= y + pParams->nDevNeighbourhood;
+ yN++) {
+ nCell = posFromXY(xN, yN, pLandscape);
+ if (nCell != _CELL_OUT_OF_RANGE) {
+ if (pLandscape->asCells[nCell].
+ nCellType == _CELL_VALID) {
+ if (pLandscape->asCells[nCell].
+ bUndeveloped == 0) {
+ testPressure++;
+ }
+ }
+ }
+ }
+ }
+ fprintf(stdout, "\t%d\n", testPressure);
+ }
+ }
+ }
+ }
+ }
+ }
+}
+
+int readParcelSizes(t_Landscape * pLandscape, t_Params * pParams)
+{
+ FILE *fIn;
+ char *szBuff;
+ int nMaxParcels;
+
+
+ pLandscape->parcelSizes = 0;
+
+ G_verbose_message("Reading patch sizes...");
+ fIn = fopen(pParams->parcelSizeFile, "rb");
+ if (fIn) {
+ szBuff = (char *)malloc(_N_MAX_DYNAMIC_BUFF_LEN * sizeof(char));
+ if (szBuff) {
+ /* just scan the file twice */
+ nMaxParcels = 0;
+ while (fgets(szBuff, _N_MAX_DYNAMIC_BUFF_LEN, fIn)) {
+ nMaxParcels++;
+ }
+ rewind(fIn);
+ if (nMaxParcels) {
+ pLandscape->aParcelSizes =
+ (int *)malloc(sizeof(int) * nMaxParcels);
+ if (pLandscape->aParcelSizes) {
+ while (fgets(szBuff, _N_MAX_DYNAMIC_BUFF_LEN, fIn)) {
+ pLandscape->aParcelSizes[pLandscape->parcelSizes] =
+ atoi(szBuff) * pParams->discountFactor;
+ if (pLandscape->aParcelSizes[pLandscape->parcelSizes]) {
+ pLandscape->parcelSizes++;
+ }
+ }
+ }
+ }
+ free(szBuff);
+ }
+ fclose(fIn);
+ }
+ return pLandscape->parcelSizes;
+}
+
+int main(int argc, char **argv)
+{
+
+ struct
+ {
+ struct Option
+ *developedFile, *devPressureFile,
+ *consWeightFile, *addVariableFiles, *nDevNeighbourhood,
+ *devpotParamsFile, *dumpFile, *outputSeries,
+ *parcelSizeFile, *discountFactor,
+ *probLookupFile,
+ *patchMean, *patchRange, *numNeighbors, *seedSearch,
+ *devPressureApproach, *alpha, *scalingFactor, *num_Regions,
+ *indexFile, *controlFileAll, *seed;
+
+ } opt;
+
+ struct
+ {
+ struct Flag *generateSeed;
+ } flg;
+
+
+ G_gisinit(argv[0]);
+
+ struct GModule *module = G_define_module();
+
+ G_add_keyword(_("raster"));
+ G_add_keyword(_("patch growing"));
+ G_add_keyword(_("urban"));
+ G_add_keyword(_("landscape"));
+ G_add_keyword(_("modeling"));
+ module->label =
+ _("Simulates landuse change using FUTure Urban-Regional Environment Simulation (FUTURES).");
+ module->description =
+ _("Module uses Patch-Growing Algorithm (PGA) to"
+ " simulate urban-rural landscape structure development.");
+
+ opt.developedFile = G_define_standard_option(G_OPT_R_INPUT);
+ opt.developedFile->key = "developed";
+ opt.developedFile->required = YES;
+ opt.developedFile->description =
+ _("Raster map of developed areas (=1), undeveloped (=0) and excluded (no data)");
+ opt.developedFile->guisection = _("Basic input");
+
+ opt.addVariableFiles = G_define_standard_option(G_OPT_R_INPUTS);
+ opt.addVariableFiles->key = "predictors";
+ opt.addVariableFiles->required = YES;
+ opt.addVariableFiles->multiple = YES;
+ opt.addVariableFiles->label = _("Names of predictor variable raster maps");
+ opt.addVariableFiles->description = _("Listed in the same order as in the development potential table");
+ opt.addVariableFiles->guisection = _("Potential");
+
+ opt.controlFileAll = G_define_standard_option(G_OPT_F_INPUT);
+ opt.controlFileAll->key = "demand";
+ opt.controlFileAll->required = YES;
+ opt.controlFileAll->description =
+ _("Control file with number of cells to convert");
+ opt.controlFileAll->guisection = _("Demand");
+
+
+ opt.devpotParamsFile = G_define_standard_option(G_OPT_F_INPUT);
+ opt.devpotParamsFile->key = "devpot_params";
+ opt.devpotParamsFile->required = YES;
+ opt.devpotParamsFile->label =
+ _("Development potential parameters for each region");
+ opt.devpotParamsFile->description =
+ _("Each line should contain region ID followed"
+ " by parameters (intercepts, development pressure, other predictors)."
+ " Values are separated by whitespace (spaces or tabs)."
+ " First line is ignored, so it can be used for header");
+ opt.devpotParamsFile->guisection = _("Potential");
+
+ opt.discountFactor = G_define_option();
+ opt.discountFactor->key = "discount_factor";
+ opt.discountFactor->type = TYPE_DOUBLE;
+ opt.discountFactor->required = YES;
+ opt.discountFactor->description = _("Discount factor of patch size");
+ opt.discountFactor->guisection = _("PGA");
+
+ opt.patchMean = G_define_option();
+ opt.patchMean->key = "compactness_mean";
+ opt.patchMean->type = TYPE_DOUBLE;
+ opt.patchMean->required = YES;
+ opt.patchMean->description =
+ _("Mean value of patch compactness to control patch shapes");
+ opt.patchMean->guisection = _("PGA");
+
+ opt.patchRange = G_define_option();
+ opt.patchRange->key = "compactness_range";
+ opt.patchRange->type = TYPE_DOUBLE;
+ opt.patchRange->required = YES;
+ opt.patchRange->description =
+ _("Range of patch compactness to control patch shapes");
+ opt.patchRange->guisection = _("PGA");
+
+ opt.numNeighbors = G_define_option();
+ opt.numNeighbors->key = "num_neighbors";
+ opt.numNeighbors->type = TYPE_INTEGER;
+ opt.numNeighbors->required = YES;
+ opt.numNeighbors->options = "4,8";
+ opt.numNeighbors->answer = "4";
+ opt.numNeighbors->description =
+ _("The number of neighbors to be used for patch generation (4 or 8)");
+ opt.numNeighbors->guisection = _("PGA");
+
+ opt.seedSearch = G_define_option();
+ opt.seedSearch->key = "seed_search";
+ opt.seedSearch->type = TYPE_INTEGER;
+ opt.seedSearch->required = YES;
+ opt.seedSearch->options = "1,2";
+ opt.seedSearch->answer="2";
+ opt.seedSearch->description =
+ _("The way location of a seed is determined (1: uniform distribution 2: development probability)");
+ opt.seedSearch->guisection = _("PGA");
+
+ opt.parcelSizeFile = G_define_standard_option(G_OPT_F_INPUT);
+ opt.parcelSizeFile->key = "patch_sizes";
+ opt.parcelSizeFile->required = YES;
+ opt.parcelSizeFile->description =
+ _("File containing list of patch sizes to use");
+ opt.parcelSizeFile->guisection = _("PGA");
+
+ opt.devPressureFile = G_define_standard_option(G_OPT_R_INPUT);
+ opt.devPressureFile->key = "development_pressure";
+ opt.devPressureFile->required = YES;
+ opt.devPressureFile->description =
+ _("Raster map of development pressure");
+ opt.devPressureFile->guisection = _("Development pressure");
+
+ opt.nDevNeighbourhood = G_define_option();
+ opt.nDevNeighbourhood->key = "n_dev_neighbourhood";
+ opt.nDevNeighbourhood->type = TYPE_INTEGER;
+ opt.nDevNeighbourhood->required = YES;
+ opt.nDevNeighbourhood->description =
+ _("Size of square used to recalculate development pressure");
+ opt.nDevNeighbourhood->guisection = _("Development pressure");
+
+ opt.devPressureApproach = G_define_option();
+ opt.devPressureApproach->key = "development_pressure_approach";
+ opt.devPressureApproach->type = TYPE_STRING;
+ opt.devPressureApproach->required = YES;
+ opt.devPressureApproach->options = "occurrence,gravity,kernel";
+ opt.devPressureApproach->description =
+ _("Approaches to derive development pressure");
+ opt.devPressureApproach->answer = "gravity";
+ opt.devPressureApproach->guisection = _("Development pressure");
+
+ opt.alpha = G_define_option();
+ opt.alpha->key = "gamma";
+ opt.alpha->type = TYPE_DOUBLE;
+ opt.alpha->required = YES;
+ opt.alpha->description =
+ _("Influence of distance between neighboring cells");
+ opt.alpha->guisection = _("Development pressure");
+
+ opt.scalingFactor = G_define_option();
+ opt.scalingFactor->key = "scaling_factor";
+ opt.scalingFactor->type = TYPE_DOUBLE;
+ opt.scalingFactor->required = YES;
+ opt.scalingFactor->description =
+ _("Scaling factor");
+ opt.scalingFactor->guisection = _("Development pressure");
+
+ opt.indexFile = G_define_standard_option(G_OPT_R_INPUT);
+ opt.indexFile->key = "subregions";
+ opt.indexFile->required = YES;
+ opt.indexFile->description = _("Raster map of subregions with categories starting with 1");
+ opt.indexFile->guisection = _("Basic input");
+
+ opt.num_Regions = G_define_option();
+ opt.num_Regions->key = "num_regions";
+ opt.num_Regions->type = TYPE_INTEGER;
+ opt.num_Regions->required = YES;
+ opt.num_Regions->description =
+ _("Number of sub-regions (e.g., counties) to be simulated");
+ opt.num_Regions->guisection = _("Basic input");
+
+ opt.probLookupFile = G_define_standard_option(G_OPT_F_INPUT);
+ opt.probLookupFile->key = "incentive_table";
+ opt.probLookupFile->required = YES;
+ opt.probLookupFile->label =
+ _("File containing incentive lookup table (infill vs. sprawl)");
+ opt.probLookupFile->description =
+ _("Format is tightly constrained. See documentation.");
+ opt.probLookupFile->guisection = _("Scenarios");
+
+ opt.consWeightFile = G_define_standard_option(G_OPT_R_INPUT);
+ opt.consWeightFile->key = "constrain_weight";
+ opt.consWeightFile->required = NO;
+ opt.consWeightFile->label =
+ _("Raster map representing development potential constraint weight for scenarios.");
+ opt.consWeightFile->description =
+ _("Values must be between 0 and 1, 1 means no constraint.");
+ opt.consWeightFile->guisection = _("Scenarios");
+
+ opt.seed = G_define_option();
+ opt.seed->key = "random_seed";
+ opt.seed->type = TYPE_INTEGER;
+ opt.seed->required = NO;
+ opt.seed->label = _("Seed for random number generator");
+ opt.seed->description =
+ _("The same seed can be used to obtain same results"
+ " or random seed can be generated by other means.");
+ opt.seed->guisection = _("Random numbers");
+
+ flg.generateSeed = G_define_flag();
+ flg.generateSeed->key = 's';
+ flg.generateSeed->label =
+ _("Generate random seed (result is non-deterministic)");
+ flg.generateSeed->description =
+ _("Automatically generates random seed for random number"
+ " generator (use when you don't want to provide the seed option)");
+ flg.generateSeed->guisection = _("Random numbers");
+
+ opt.dumpFile = G_define_standard_option(G_OPT_R_OUTPUT);
+ opt.dumpFile->key = "output";
+ opt.dumpFile->required = YES;
+ opt.dumpFile->description =
+ _("State of the development at the end of simulation");
+ opt.dumpFile->guisection = _("Output");
+
+ opt.outputSeries = G_define_standard_option(G_OPT_R_BASENAME_OUTPUT);
+ opt.outputSeries->key = "output_series";
+ opt.outputSeries->required = NO;
+ opt.outputSeries->label =
+ _("Basename for raster maps of development generated after each step");
+ opt.outputSeries->guisection = _("Output");
+ // TODO: add mutually exclusive?
+ // TODO: add flags or options to control values in series and final rasters
+
+
+ // provided XOR generated
+ G_option_exclusive(opt.seed, flg.generateSeed, NULL);
+ G_option_required(opt.seed, flg.generateSeed, NULL);
+
+ if (G_parser(argc, argv))
+ exit(EXIT_FAILURE);
+
+ long seed_value;
+
+ if (flg.generateSeed->answer) {
+ seed_value = G_srand48_auto();
+ G_message("Generated random seed (-s): %ld", seed_value);
+ }
+ if (opt.seed->answer) {
+ seed_value = atol(opt.seed->answer);
+ // this does nothing since we are not using GRASS random function
+ G_srand48(seed_value);
+ G_message("Read random seed from %s option: %ld",
+ opt.seed->key, seed_value);
+ }
+ // although GRASS random function is initialized
+ // the general one must be done separately
+ // TODO: replace all by GRASS random number generator?
+ srand(seed_value);
+
+ // TODO: move this back to local variables
+ //t_Params sParams;
+ //t_Landscape sLandscape;
+
+ /* blank everything out */
+ memset(&sParams, 0, sizeof(t_Params));
+
+ /* parameters dependednt on region */
+ sParams.xSize = Rast_window_rows();
+ sParams.ySize = Rast_window_cols();
+ G_message(_("Running on %d rows and %d columns (%d cells)"),
+ sParams.xSize, sParams.ySize, sParams.xSize * sParams.ySize);
+
+ /* set up parameters */
+ sParams.developedFile = opt.developedFile->answer;
+ sParams.devPressureFile = opt.devPressureFile->answer;
+ sParams.consWeightFile = opt.consWeightFile->answer;
+ sParams.numAddVariables = 0;
+
+ size_t num_answers = 0;
+ while (opt.addVariableFiles->answers[num_answers]) {
+ sParams.addVariableFile[num_answers] = opt.addVariableFiles->answers[num_answers];
+ num_answers++;
+ }
+ sParams.numAddVariables = num_answers;
+ // TODO: dyn allocate file list
+ sParams.nDevNeighbourhood = atoi(opt.nDevNeighbourhood->answer);
+
+ sParams.dumpFile = opt.dumpFile->answer;
+ sParams.outputSeries = opt.outputSeries->answer;
+
+ sParams.parcelSizeFile = opt.parcelSizeFile->answer;
+
+ sParams.discountFactor = atof(opt.discountFactor->answer);
+
+ // TODO: remove all sortProbs != 0 code if it does not make sense for Stochastic 2
+ // always 0 for Stochastic 2
+ sParams.sortProbs = 0;
+
+ // TODO: remove remaining stochastic 1 and deterministic code
+ // stochastic 1 and deterministic not maintained, reimplementing considered as easier
+ // then fixing it (old code will be still available in the old version)
+ // always use Stochastic II
+ sParams.nAlgorithm = _N_ALGORITHM_STOCHASTIC_II;
+
+ if (sParams.nAlgorithm == _N_ALGORITHM_STOCHASTIC_II) {
+
+ int parsedOK, i;
+ FILE *fp;
+ char inBuff[N_MAXREADINLEN];
+ char *pPtr;
+
+ G_verbose_message("Reading probability lookup ...");
+ sParams.probLookupFile = opt.probLookupFile->answer;
+
+ fp = fopen(sParams.probLookupFile, "r");
+ if (fp) {
+ parsedOK = 0;
+ if (fgets(inBuff, N_MAXREADINLEN, fp)) {
+ if (inBuff[0] == ',') {
+ sParams.nProbLookup = atoi(inBuff + 1);
+ if (sParams.nProbLookup > 0) {
+ sParams.adProbLookup =
+ (double *)malloc(sizeof(double) *
+ sParams.nProbLookup);
+ if (sParams.adProbLookup) {
+ parsedOK = 1;
+ i = 0;
+ while (parsedOK && i < sParams.nProbLookup) {
+ parsedOK = 0;
+ if (fgets(inBuff, N_MAXREADINLEN, fp)) {
+ if (pPtr = strchr(inBuff, ',')) {
+ parsedOK = 1;
+ sParams.adProbLookup[i] =
+ atof(pPtr + 1);
+ G_debug(3,
+ "probability lookup table: i=%d, i/(n-1)=%f, result=%f",
+ i,
+ i * 1.0 /
+ (sParams.nProbLookup - 1),
+ sParams.adProbLookup[i]);
+ }
+ }
+ i++;
+ }
+ }
+ }
+ }
+ }
+ if (!parsedOK) {
+ G_fatal_error("Error parsing probability lookup file '%s'",
+ sParams.probLookupFile);
+ }
+ fclose(fp);
+ }
+ else {
+ G_fatal_error("Error opening probability lookup file '%s'",
+ sParams.probLookupFile);
+ }
+
+ sParams.patchMean = atof(opt.patchMean->answer);
+ sParams.patchRange = atof(opt.patchRange->answer);
+ sParams.numNeighbors = atoi(opt.numNeighbors->answer);
+ // TODO: convert to options or flag: 1: uniform distribution 2: based on dev. proba.
+ sParams.seedSearch = atoi(opt.seedSearch->answer);
+ sParams.devPressureApproach = atoi(opt.devPressureApproach->answer);
+ if (strcmp(opt.devPressureApproach->answer, "occurrence") == 0)
+ sParams.devPressureApproach = 1;
+ else if (strcmp(opt.devPressureApproach->answer, "gravity") == 0)
+ sParams.devPressureApproach = 2;
+ else if (strcmp(opt.devPressureApproach->answer, "kernel") == 0)
+ sParams.devPressureApproach = 3;
+ else
+ G_fatal_error(_("Approach doesn't exist"));
+ if (sParams.devPressureApproach != 1) {
+ sParams.alpha = atof(opt.alpha->answer);
+ sParams.scalingFactor = atof(opt.scalingFactor->answer);
+ }
+
+ sParams.num_Regions = atoi(opt.num_Regions->answer);
+ sParams.indexFile = opt.indexFile->answer;
+ sParams.controlFileAll = opt.controlFileAll->answer;
+ }
+
+ if (sParams.num_Regions > 1)
+ readDevPotParams(&sParams, opt.devpotParamsFile->answer);
+
+ readDevDemand(&sParams);
+ /* allocate memory */
+ if (buildLandscape(&sLandscape, &sParams)) {
+ /* read data */
+ if (readData(&sLandscape, &sParams)) {
+ readData4AdditionalVariables(&sLandscape, &sParams);
+ readIndexData(&sLandscape, &sParams);
+ if (readParcelSizes(&sLandscape, &sParams)) {
+ //testDevPressure(&sLandscape, &sParams);
+ /* do calculation and dump result */
+ //one region, enable the following
+ //updateMap(&sLandscape, &sParams);
+
+ //multiple regions
+ updateMapAll(&sLandscape, &sParams);
+ }
+ else {
+ G_fatal_error("Reading patch sizes failed");
+ }
+ }
+ else {
+ G_fatal_error("Reading input maps failed");
+ }
+ /* could put in routines to free memory, but OS will garbage collect anyway */
+ }
+ else {
+ G_fatal_error("Initialization failed");
+ }
+
+ return EXIT_SUCCESS;
+}
+
+int getUnDevIndex(t_Landscape * pLandscape)
+{
+ float p = rand() / (double)RAND_MAX;
+ int i;
+
+ for (i = 0; i < pLandscape->undevSites; i++) {
+ if (p < pLandscape->asUndev[i].cumulProb) {
+ return i;
+ }
+ }
+ // TODO: returning at least something but should be something more meaningful
+ return 0;
+}
+
+int getUnDevIndex1(t_Landscape * pLandscape, int regionID)
+{
+ float p = rand() / (double)RAND_MAX;
+
+ G_debug(3, _("getUnDevIndex1: regionID=%d, num_undevSites=%d, p=%f"),
+ regionID, pLandscape->num_undevSites[regionID], p);
+ int i;
+
+ for (i = 0; i < pLandscape->num_undevSites[regionID]; i++) {
+ if (p < pLandscape->asUndevs[regionID][i].cumulProb) {
+ G_debug(3, _("getUnDevIndex1: cumulProb=%f"),
+ pLandscape->asUndevs[regionID][i].cumulProb);
+ return i;
+ }
+ }
+ // TODO: returning at least something but should be something more meaningful
+ return 0;
+}
+
+void findAndSortProbsAll(t_Landscape * pLandscape, t_Params * pParams,
+ int step)
+{
+ int i, lookupPos;
+ t_Cell *pThis;
+ int id;
+
+ /* update calcs */
+ G_verbose_message("Recalculating probabilities");
+ for (i = 0; i < pParams->num_Regions; i++) {
+ pLandscape->num_undevSites[i] = 0;
+ }
+ float val = 0.0;
+
+ for (i = 0; i < pLandscape->totalCells; i++) {
+ pThis = &(pLandscape->asCells[i]);
+ if (pThis->nCellType == _CELL_VALID) {
+ if (pThis->bUndeveloped) {
+ if (pThis->consWeight > 0.0) {
+ id = pThis->index_region - 1;
+ if (pThis->index_region == -9999)
+ continue;
+ /* note that are no longer just storing the logit value,
+ but instead the probability
+ (allows consWeight to affect sort order) */
+
+ if (id < 0 || id >= pParams->num_Regions)
+ G_fatal_error(_("Index of region %d is out of range [0, %d] (region is %d)"),
+ id, pParams->num_Regions - 1,
+ pThis->index_region);
+
+ if (pLandscape->num_undevSites[id] >=
+ pLandscape->asUndevs[id].size())
+ pLandscape->asUndevs[id].resize(pLandscape->
+ asUndevs[id].size() *
+ 2);
+
+ pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].
+ cellID = i;
+ val = getDevProbability(pThis, pParams);
+ pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].
+ logitVal = val;
+ G_debug(2, "logit value %f", val);
+ pThis->devProba = val;
+ /* lookup table of probabilities is applied before consWeight */
+ if (pParams->nAlgorithm == _N_ALGORITHM_STOCHASTIC_II) {
+ /* replace with value from lookup table */
+ lookupPos = (int)(pLandscape->asUndevs[id]
+ [pLandscape->num_undevSites[id]].
+ logitVal * (pParams->nProbLookup -
+ 1));
+ if (lookupPos >= pParams->nProbLookup ||
+ lookupPos < 0)
+ G_fatal_error
+ ("lookup position (%d) out of range [0, %d]",
+ lookupPos, pParams->nProbLookup - 1);
+ pLandscape->asUndevs[id][pLandscape->
+ num_undevSites[id]].
+ logitVal = pParams->adProbLookup[lookupPos];
+ }
+ pThis->devProba =
+ pLandscape->asUndevs[id][pLandscape->
+ num_undevSites[id]].logitVal;
+ // discount by a conservation factor
+ pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].logitVal *= pThis->consWeight;
+ /* need to store this to put correct elements near top of list */
+ pLandscape->asUndevs[id][pLandscape->num_undevSites[id]].bUntouched = pThis->bUntouched;
+ if (pLandscape->asUndevs[id]
+ [pLandscape->num_undevSites[id]].logitVal > 0.0) {
+ /* only add one more to the list if have a positive probability */
+ pLandscape->num_undevSites[id]++;
+ }
+ }
+ }
+ }
+ }
+ /* sort */
+ /* can only be zero for algorithm=stochastic_ii */
+ if (pParams->sortProbs) {
+ G_debug(1, "Sorting %d unconserved undeveloped sites",
+ pLandscape->undevSites);
+ qsort(pLandscape->asUndev, pLandscape->undevSites, sizeof(t_Undev),
+ undevCmpReverse);
+ }
+ else {
+ G_debug(1, "Skipping sort as choosing cells randomly");
+ }
+ //calculate cumulative probability // From Wenwu Tang
+ int j;
+
+ for (j = 0; j < pParams->num_Regions; j++) {
+ double sum = pLandscape->asUndevs[j][0].logitVal;
+
+ for (i = 1; i < pLandscape->num_undevSites[j]; i++) {
+ pLandscape->asUndevs[j][i].cumulProb =
+ pLandscape->asUndevs[j][i - 1].cumulProb +
+ pLandscape->asUndevs[j][i].logitVal;
+ sum = sum + pLandscape->asUndevs[j][i].logitVal;
+ }
+ for (i = 0; i < pLandscape->num_undevSites[j]; i++) {
+ pLandscape->asUndevs[j][i].cumulProb =
+ pLandscape->asUndevs[j][i].cumulProb / sum;
+ }
+ }
+}
Copied: grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html (from rev 65540, grass-addons/grass7/raster/r.futures/r.futures.html)
===================================================================
--- grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html (rev 0)
+++ grass-addons/grass7/raster/r.futures/r.futures.pga/r.futures.pga.html 2015-07-06 18:56:57 UTC (rev 65541)
@@ -0,0 +1,157 @@
+<h2>DESCRIPTION</h2>
+
+Module <em>r.futures.pga</em> is part of <a href="r.futures.html">FUTURES</a>
+land change model.
+This module uses stochastic Patch-Growing Algorithm (PGA)
+and a combination of field-based and object-based representations
+to simulate land changes.
+
+PGA simulates undeveloped to developed land change by iterative site selection
+and a contextually aware region growing mechanism.
+Simulations of change at each time step feed development pressure back
+to the POTENTIAL submodel, influencing site suitability for the next step.
+
+<h3>Patch growing</h3>
+Patches are constructed in
+three steps. First, a potential seed is randomly selected from available cells.
+In case <b>seed_search</b> is 2, the probability value (based on POTENTIAL)
+of the seed is tested using Monte Carlo approach, and if it doesn't survive,
+new potential seed is selected and tested.
+Second, using a 4- or 8-neighbor (see <b>num_neighbors</b>) search rule PGA grows the patch.
+PGA decides on the suitability of contiguous cells based on their
+underlying development potential and distance to the seed adjusted
+by compactness parameter given in <b>compactness_mean</b> and <b>compactness_range</b>.
+The size of the patch is determined by randomly selecting a patch size from <b>patch sizes</b> file
+and multiplied by <b>discount_factor</b>. To find optimal values
+for patch sizes and compactness, use module <a href="r.futures.calib.html">r.futures.calib</a>.
+Once a cell is converted, it remains developed.
+PGA continues to grow patches until the per capita land demand is satisfied.
+
+<h3>Development pressure</h3>
+Development pressure is a dynamic spatial variable
+derived from the patch-building process of PGA and associated with the POTENTIAL submodel.
+At each time step, PGA updates the POTENTIAL probability surface based on land change,
+and the new development pressure then affects future land change.
+The initial development pressure is computed using module <a href="r.futures.devpressure.html">
+r.futures.devpressure</a>. The same input parameters of this module
+(<b>gamma</b>, <b>scaling factor</b> and <b>n_dev_neighbourhood</b>)
+are then used as input for <a href="r.futures.pga.html">
+r.futures.pga</a>.
+
+<h3>Scenarios</h3>
+Scenarios involving policies that encourage infill versus sprawl
+can be explored using the <b>incentive_table</b> parameter,
+which uses a power function to transform the probability in POTENTIAL.
+
+A simple incentive table can be generated using the
+following commands in Python:
+<div class="code"><pre>
+import numpy as np
+power = 2
+x = np.linspace(0, 1, 1001)
+np.savetxt(fname="incentive_table.txt", fmt='%1.3f', delimiter=',',
+ X=np.column_stack((x, np.power(x, power))))
+</pre></div>
+
+<center>
+<img src="incentive.png">
+<p>
+Figure: Effect of incentive table on development probability
+</center>
+<p>
+Additionally, parameter <b>constrain_weight</b> (raster map from 0 to 1)
+enables to include policies (such as new regulations or fees) which limit the
+development in certain areas.
+The probability surface is simply multiplied by the <b>constrain_weight</b> values,
+which results in decreased site suitability in areas,
+where the <b>constrain_weight</b> values are lower than 1.
+
+<h3>Output</h3>
+After the simulation ends, raster specified in parameter <b>output</b> is written.
+If optional parameter <b>output_series</b> is specified, additional output
+is a series of raster maps for each step.
+Cells with value 0 represents the initial development, values >= 1 then represent
+the step in which the cell was developed. Undeveloped cells have value -1.
+<p>
+<center>
+<img src="r_futures.png">
+<p>
+Figure: Output map of developed areas
+</center>
+
+<p>
+<center>
+<img src="r_futures_detail.png">
+<p>
+Figure: Detail of output map
+</center>
+
+
+<h2>EXAMPLE</h2>
+
+<div class="code"><pre>
+r.futures.pga -s developed=lc96 predictors=d2urbkm,d2intkm,d2rdskm,slope demand=demand.txt devpot_params=devpotParams.csv discount_factor=0.6 compactness_mean=0.4 compactness_range=0.08 num_neighbors=4 seed_search=2 patch_sizes=patch_sizes.txt development_pressure=gdp n_dev_neighbourhood=10 development_pressure_approach=gravity gamma=2 scaling_factor=1 subregions=subregions num_regions=9 incentive_table=probLookup.csv constrain_weight=weight_1 output=final_results output_series=development
+</pre></div>
+
+
+<h2>KNOWN ISSUES</h2>
+
+Module holds all data in the memory, so for large areas the memory requirements
+are high.
+
+Some inputs are not yet documented.
+
+
+<h2>REFERENCES</h2>
+
+<ul>
+<li>
+ Meentemeyer, R. K., Tang, W., Dorning, M. A., Vogler, J. B., Cunniffe, N. J., & Shoemaker, D. A. (2013).
+ <a href="http://dx.doi.org/10.1080/00045608.2012.707591">FUTURES: Multilevel Simulations of Emerging
+ Urban-Rural Landscape Structure Using a Stochastic Patch-Growing Algorithm</a>.
+ Annals of the Association of American Geographers, 103(4), 785-807.
+<li>Dorning, M. A., Koch, J., Shoemaker, D. A., & Meentemeyer, R. K. (2015).
+ <a href="http://dx.doi.org/10.1016/j.landurbplan.2014.11.011">Simulating urbanization scenarios reveals
+ tradeoffs between conservation planning strategies</a>.
+ Landscape and Urban Planning, 136, 28-39.</li>
+</ul>
+
+
+<h2>SEE ALSO</h2>
+
+<a href="r.futures.html">FUTURES</a>,
+<em><a href="r.futures.devpressure.html">r.futures.devpressure</a></em>,
+<em><a href="r.futures.calib.html">r.futures.calib</a></em>,
+<em><a href="r.sample.category.html">r.sample.category</a></em>
+
+
+<h2>AUTHORS</h2>
+
+<p>
+<em>Corresponding author:</em><br>
+Ross K. Meentemeyer, rkmeente ncsu edu,
+<a href="http://geospatial.ncsu.edu/">Center for Geospatial Analytics, NCSU</a>
+
+<p>
+<em>Original standalone version:</em><br>
+
+Ross K. Meentemeyer *<br>
+Wenwu Tang *<br>
+Monica A. Dorning *<br>
+John B. Vogler *<br>
+Nik J. Cunniffe *<br>
+Douglas A. Shoemaker *<br>
+Jennifer A. Koch **<br>
+
+<br>
+* Department of Geography and Earth Sciences, UNC Charlotte<br>
+** <a href="http://geospatial.ncsu.edu/">Center for Geospatial Analytics, NCSU</a><br>
+
+<p>
+<em>Port to GRASS GIS and GRASS-specific additions:</em><br>
+
+Vaclav Petras, <a href="http://geospatial.ncsu.edu/osgeorel/">NCSU OSGeoREL</a><br>
+
+
+<p>
+<i>Last changed: $Date$</i>
Copied: grass-addons/grass7/raster/r.futures/r.futures.pga/r_futures.png (from rev 65540, grass-addons/grass7/raster/r.futures/r_futures.png)
===================================================================
(Binary files differ)
Copied: grass-addons/grass7/raster/r.futures/r.futures.pga/r_futures_detail.png (from rev 65540, grass-addons/grass7/raster/r.futures/r_futures_detail.png)
===================================================================
(Binary files differ)
Deleted: grass-addons/grass7/raster/r.futures/r_futures.png
===================================================================
(Binary files differ)
Deleted: grass-addons/grass7/raster/r.futures/r_futures_detail.png
===================================================================
(Binary files differ)
More information about the grass-commit
mailing list