UNIVERSIDAD COMPLUTENSE DE MADRID FACULTAD DE GEOGRAFIA E HISTORIA MASTER UNIVERSITARIO EN TECNOLOGIAS DE LA INFORMACION GEOGRAFICA TRABAJO FIN DE MASTER CURSO: 2020-2021 Accesibilidad dinámica en modelos de Autómata Celular (AC) de usos del suelo urbano. Dynamic accesibility on Cellular Automata (CA) urban land use models. PEREZ-BELLO PIÑEIRO, JULIO-JOSE CONVOCATORIA: SEPTIEMBRE DIRECTORES: Richard James Hewitt & Jaime Salvador Díaz Pacheco. Departamento de Geografía e Historia. AGRADECIMIENTOS A Richard James Hewitt y a Jaime Salvador Díaz Pacheco, por enseñarme el 18 de diciembre de 2020 que es un Autómata Celular y sus posibles aplicaciones en geografía. Un simple agradecimiento hacia ti, Richard, me parece muy egoísta por mi parte. Has sido la persona de la que más he aprendido y más horizontes (que desconocida completamente) me ha abierto con la introducción en el lenguaje de programación R utilizado como un GIS. Hace un año no me hubiera creído que acabaría realizando un trabajo de GIS en un lenguaje de programación. Esto no es solo un trabajo final de máster, es una enseñanza paso a pasco desde 0, sin tener idea de programar en R, por el cual hubo un proceso de aprendizaje por mi parte por introducirme en otro sistema completamente diferenciado y nuevo para mi sin pasar por ArcGIS o QGIS. Es un logro personal. A mis padres por haberme dado la oportunidad de seguir formándome y confiado en mí todo este tiempo y apoyado desde Galicia, brindándome todo lo que han podido y más. Poco más puedo decir de unos padres que me lo han dado todo, sin tener que dar explicaciones. Simplemente: Hazlo y aprovecha tu oportunidad mientras se pueda, que muchas de ellas solo ocurren una vez. Los años pasan y nuestras consecuencias tienen una acción. Intenta por lo menos sacar lo positivo de todo aquello que hagas, sea bueno o no. Algo siempre se aprende. RESUMEN Los Autómatas Celulares aplicados a los sistemas de información geográfica han permitido el uso de nuevos paradigmas científicos en el análisis del medio. El presente trabajo se emplaza dentro de la geografía urbana, donde a partir de una serie de reglas aplicadas sobre los Autómatas Celulares, como son la vecindad, la accesibilidad, el relieve, el factor estocástico y la aptitud, se iniciará el análisis del medio urbano de la Comunidad de Madrid. El análisis se realizará a partir de dos tipos de modelos: estático y dinámico, siendo el lenguaje de programación R el método por el cual se crearán los Scripts y las simulaciones finales. A partir de una serie de calibraciones iniciales, se jugará con el Autómata con el fin de conseguir unos resultados óptimos y lo más reales posibles. Tras ello, se comenzará a usar los métodos anteriormente dichos y se compararán con diferentes herramientas como MCK o Fragstats. El estudio pretende analizar como la accesibilidad ( viaria en este caso), junto con una vecindad previamente establecida por el propio autor, conocer como el uso del suelo urbano se desarrolla desde el año 2000 hasta el año 2050 sobre la Comunidad de Madrid con la aplicación de una accesibilidad dinámica y estática. El objetivo final, por tanto, es observar la evolución del suelo urbano en función de las variables de accesibilidad y vecindad, o como el desarrollo del suelo sobre la CAM genera nuevas estructuras urbanas a partir de unos datos iniciales de suelo a partir de la capa de Madrid Land Use. Los resultados serán positivos en este aspecto sobre los modelos analizados, pudiéndose cerciorar que ambos modelos para el caso de estudio son completamente válidos. Como implicaciones futuras, se puede llegar a generar una gran potencialidad sobre este Script y los propios Autómatas a la hora de crear nuevas simulaciones, no solo de suelo urbanos, sino de cualquier paradigma aplicado, como por ejemplo las migraciones de aves o las áreas de distribución de diferentes especies animales a partir de un patrón más simple. Palabras Clave: R, MCK, SIMLANDER, Madrid Land Use, Autómatas Celulares. ABSTRACT Cellular Automata applied to geographic information systems have allowed the use of new scientific paradigms in the analysis of the environment. The present work is located within urban geography, where, based on a series of rules applied to Cellular Automata, such as neighborhood, accessibility, relief, stochastic factor and fitness, the analysis of the urban environment will begin. of the Community of Madrid. The analysis will be carried out from two types of models: static and dynamic, being the R programming language the method by which the final scripts and simulations will be created. Starting from a series of initial calibrations, the Automaton will be played in order to achieve optimal results and as realistic as possible. After that, the aforementioned methods will begin to be used and they will be compared with different tools such as MCK or Fragstats. The study aims to analyze how accessibility (road in this case), together with a neighborhood previously established by the author himself, to know how the use of urban land develops from 2000 to 2050 on the Community of Madrid with the application of a dynamic and static accessibility. The final objective, therefore, is to observe the evolution of urban land as a function of the variables of accessibility and neighborhood, or how the development of the land on the CAM generates new urban structures based on initial soil data from the layer of Madrid Land Use. The results will be positive in this regard on the models analyzed, being able to make sure that both models for the case study are completely valid. As future implications, a great potential can be generated on this Script and the Automata themselves when creating new simulations, not only of urban land, but of any applied paradigm, such as migrations of birds or areas of distribution of different animal species from a simpler pattern. Keywords: R, MCK, SIMLANDER, Madrid Land Use, Cellular Automata. 1 INDICE DE CONTENIDOS 1 – INTRODUCCION .............................................................................. 2 2 - MARCO TEORICO ........................................................................... 3 2.1-AUTÓMATA CELULAR (AC) Y GEOGRAFÍA. ........................................ 3 2.1.1 - El juego de la vida; por John Conway. ...................................... 5 2.1.2 - La importancia de los fractales. ................................................. 6 3 – METODOLOGIA .............................................................................. 8 3.1 MÉTODOS DE APLICACIÓN DEL AUTÓMATA CELULAR. ......................... 8 3.1.1 – Accesibilidad ............................................................................ 8 3.1.2 - Factor aleatorio ....................................................................... 10 3.1.2.1 - Establecer los diferentes tipos de uso del suelo sobre el área de estudio: ............................................................................................................................. 10 3.1.2.2 - La escala: ............................................................................................... 10 3.1.3 - Vecindad ................................................................................. 12 3.1.3.1 – Representación de los usos urbanos en el modelo:............................... 12 3.1.4 – Aptitud ................................................................................... 14 3.1.5 - Cálculo del potencial de transición .......................................... 17 4 - SIMLANDER .................................................................................... 18 4.1 - EL USO DE SOFTWARE R COMO GIS. ............................................... 19 4.2-APLICACIÓN DE SIMLANDER MEDIANTE R (CALIBRACIÓN DEL SISTEMA)................................................................................................ 19 4.2.1 Accesibilidad estática (2000-2009). ........................................... 21 4.2.1.1 - Comparativa de los mapas estáticos entre los años 2000-2009. ........... 26 4.3 - MODELIZACIÓN DINÁMICA Y ESTÁTICA ENTRE LOS AÑOS 2009-2050. .............. 30 4.3.1 - Comparación con Map Comparison Kit. ............................................ 36 4.3.2 - Comparación con Fragstats. ............................................................ 40 5 - CONCLUSIONES ............................................................................. 42 6 – REFERENCIAS BIBLIOGRAFICAS ............................................ 43 7-ENLACES WEB ................................................................................. 45 8- ANEXOS ............................................................................................ 46 2 1 – INTRODUCCION La Geografía es una ciencia social, de la cual debemos tener en cuenta 3 ideas principales. La plasmación de la realidad física, humana, regional, no solo a través del estudio gráfico o sintáctico, sino también analíticamente. Las infinitas posibilidades que han emergido desde principios de los 90 con la liberación (Open Source) de las Tecnologías de la Información Geográfica. Las TIG han generado en el ámbito geográfico un nuevo sinfín de posibilidades de dispersión en cuestiones científicas, aplicando la computación. Estos hechos, han propiciado la consecución de crear nuevos horizontes científicos, como es el caso de los estudios dedicados a los usos del suelo aplicados a los formatos ráster y vectorial. Estos dos tipos de formatos han conseguido recrear las situaciones geográficas digitalmente y en dos ámbitos simples, pero a la vez diferenciados. Las herramientas SIG han aportado un nuevo mundo digital, destacando en nuestro asunto, la metamorfosis de las ciudades a través de cálculos matemáticos aplicados sobre una serie de reglas algebraicas a los usos (Tobler 1979 & Coucelis 1980). Las posibilidades de concebir a través de una serie de procesos matemáticos la generación de una ciudad desde su nacimiento a través de un píxel, nos permite averiguar sus defectos, su crecimiento regular o irregular, o por qué la ciudad crece más por el este que por el oeste. Todos estos píxeles que se van autogenerando a través de las reglas de juego establecidas, son lo que denominamos Autómatas Celulares. En un SIG, los Autómatas Celulares son entes en formato ráster (píxeles) que se van generando en función de la probabilidad estimada para cada celda o píxel, según el Autómata va desarrollándose. A medida que se desarrolla, el Autómata va calculando en cada simulación o generación, esa probabilidad de manera autónoma para cada célula raster. Su principal virtud radica en la interacción entre las células de manera sencilla, manifestando un comportamiento múltiple globalmente, debido a que a partir de unas reglas relativamente simples se generan estructuras complejas. Sus propiedades comienzan siendo muy simples desde su creación como píxel individual, sin embargo, a la vez que este píxel nace, su alrededor va calculando nuevas posibilidades de crecimiento sobre el uso del suelo y, por tanto, determinando como serán los píxeles que circundarán este primero. Nos permiten simular la evolución del espacio geográfico. Así comienza el desarrollo de este tipo de autómatas, con ciertos rasgos aplicables a la teoría celular: Nacimiento, creación, destrucción. Nuestra área de estudio se basará en conocer como la comunidad de Madrid se adecúa a las posibilidades de cambios a través de las redes viarias. Los modelos Autómata Celular de usos del suelo urbano (White et al. SLEUTH, ArcModiland, SIMLANDER) habitualmente incorporan la accesibilidad, normalmente las redes viarias o ferroviarias, como un factor determinante en las simulaciones de la evolución espacio-temporal de las ciudades. 3 Está bien documentado en la literatura científica el papel clave que juega la movilidad en el desarrollo de la estructura espacial de las zonas urbanas, de tal manera que las infraestructuras pueden conducir una ciudad hacia una forma más compacta, o bien más dispersa, en función de su ubicación y su tipo (vehículos de motor o bicicletas, transporte público o privado) y su tamaño (autovía o carretera local). Sin embargo, en los modelos de usos del suelo urbano no es habitual tener en cuenta el dinamismo de las propias redes, y las simulaciones típicamente se desarrollan en función de un mapa base de accesibilidad, derivado del trazado de la red de infraestructuras al principio de la simulación. Hay pocos los trabajos que tratan de simular la evolución de la red de transporte en función de la forma urbana, y de la forma urbana en función de la red de transporte, a la misma vez. En una simulación de este tipo las infraestructuras cambian su ubicación, tipo o tamaño según el patrón simulado, influyendo a su vez, el patrón urbano en el siguiente paso temporal. El presente trabajo pretende indagar sobre este tema, a través de una serie de experimentos llevados a cabo con la base de datos Madrid Land Use (MLU), y los entornos de simulación urbana de AC, ArcModiland (Díaz Pacheco & Moya Gómez 2013) y SIMLANDER (Hewitt et al. 2013). 2 - MARCO TEORICO 2.1-Autómata Celular (AC) y Geografía. Cuando hablamos de AC, debemos atenernos al uso del formato ráster en Geografía. Usamos este tipo de formato con este término por el interés en observar los fenómenos de cambio que se producen sobre los usos del suelo. Los formatos ráster, son originalmente en su estructura más simples, elementos en formato de píxeles de representación numérica, alfanumérica o claves algebraicas en formato cartográfico. Con este tipo de formato, se pueden elaborar simulaciones basadas en una serie propiedades previamente definidas, a través del uso de la geoestadística, la probabilidad o el uso de algoritmos matemáticos previamente generados mediante Inteligencia Artificial (IA) para recrear su procreación (Bhatta 2009). Su campo de aplicación es infinito (se ha declarado como “A Turing machine”, es decir, un sistema computacional universal de aplicación por el hecho de que cualquier problema que se puede expresar en forma de algoritmo, el autómata puede resolverlo) incluso en el campo de la teoría celular basado en su reproducción y movimiento aleatorio, siendo un punto de inflexión clave y de unión entre la biología y la geografía basada principalmente en los procesos geoestadísticos. Es aquí cuando debemos tener en cuenta siempre la Ley de Tobler en la cual, cualquier elemento que se encuentre cercano, notará los cambios de una manera más similar y acorde, que aquellos que se encuentren más alejados, siendo este cambio, progresivo con el elemento espacio-distancia. Nuestro caso de aplicación de la teoría celular se basa en como las ciudades, formulados como entes de vida, van creándose, modificándose, 4 desapareciendo o incluso renaciendo o uniéndose a sus semejantes por su cercanía (Portugali et al. 2021). En el campo de la Geografía sería su vecindad más próxima la que dictamine las reglas de crecimiento en ciertas aplicaciones basadas en medioambiente y el área de urbanismo. No nos basamos solo en la ciudad como un ciclo de vida, sino también de atracción y de influencias sobre aquellos entes similares o superiores que de una forma u otra generan cambios espaciales generales e incluso reproducen la vida. Estos entes que parecen poseer una vida, los llamaremos autómatas celulares. A pesar de que no sea posiblemente un vocablo demasiado extendido para la mayoría de la población, en el campo de las ciencias es conocido, siendo desarrollado a finales de los años 40 por los físicos John Von Neumann y Stanislaw Ulam. Su principal característica es la extrema simpleza que posee la célula en un estadio inicial. Conforme se altera y transforma al crecer, la disposición de sus átomos (píxeles) va difiriendo de la inicial y es aquí donde radica la complejidad de la célula autómata: la aleatoriedad constante sobre cómo se desarrollará. Las primeras aplicaciones en el campo del urbanismo se remontan a los investigadores Cechini y Viola (1993), quienes lo adaptaron a los modelos urbanos de crecimiento. Vieron como las ciudades nacen, se reproducen y crecen en conjunto sobre un azar constante (perturbaciones estocásticas), el cual puede llegar a ser calculado utilizando un gran número de parámetros a tener en cuenta, que, hasta cierto punto, imposibilitan una implementación correcta, puesto que se necesita programarlas de una manera exacta para cada tipo de pixel. Es por ello por lo que se trata de un proceso lógicamente estocástico propio del azar en algunos casos y por el cual se usa un factor aleatorio: para dotar al sistema de una variación en las condiciones iniciales que generen cambios y diferenciaciones en los experimentos (Lorenz 1962). Profundizando en los procesos de uso del suelo sobre un modelo celular debemos tener una serie de inflexiones como el tipo de tamaño del ráster, el radio de vecindad que utilizará, en cualquier caso, pues determinará en otras circunstancias un desarrollo diferente sobre ello. Si realizáramos una simulación a partir de unos parámetros predefinidos sobre un área de una ciudad de grandes dimensiones, veríamos como los centros celulares seguirían un patrón homogéneo siendo totalmente opuesto a los lugares más alejados del centro. Esto es consecuencia de la interacción con otras células, en este caso ráster o píxeles de diferente influencia sobre ellas, que generan procesos aleatorios en el extrarradio de una ciudad en crecimiento, en función del cálculo de la vecindad usado. Por tanto, el factor principal y más importante, es el estadio inicial de un ente en el cual comenzara a redirigir su crecimiento por vecindad ráster, siguiendo una serie de patrones espacio-temporales (Lorenz 1962). 5 2.1.1 - El juego de la vida; por John Conway. Uno de los conceptos más importantes a conocer sobre el funcionamiento de los Autómatas celulares, independientemente de si su creación es a partir de las reglas del juego que establezcamos, es que debemos tener en cuenta cómo funciona internamente a semejanza de un ente vivo, siendo el juego de la vida un buen ejemplo didáctico para ello. Como hemos dicho anteriormente, los Autómatas se comportan como entidades vivas, por ende, se transforman y se mueven en el espacio y tiempo a través de una serie de reglas matemáticas, descritas por el matemático británico John Conway en 1970. A pesar de la complejidad de los patrones que se pueden generar, las reglas son muy simples. En su caso, establece que la configuración primaria debe atenerse a que el crecimiento no sea limitado sino variable, así como cambios durante su crecimiento generando alteraciones profundas sobre la configuración a medida que se va desarrollando. El principal reto, es que el comportamiento de este crecimiento sea, no predecible, evitando por un lado que genere algo predeterminado. Por ello, el patrón inicial que cada píxel o célula inicial presenta es al igual que en el caso de la aplicación en geografía y la teoría de Tobler, que aquellas células adyacentes afectan sobre el estadio inicial, de tal manera que las celdas contiguas impidan una regularidad de crecimiento. Así se generan áreas con crecimientos, eliminaciones o nacimientos celulares. Conway nos indica que el patrón no puede crecer sin límite, sino que la configuración inicial no puede tener un numero infinito de movimientos a la hora de crearse. El único modo de ser un sistema infinito sería generando nuevos cálculos sucesivos, pero, debemos tener en cuenta que a nivel celular si el estadio inicial es de 1 a 3 células, y estas no poseen una vecindad próxima, irremediablemente perecerán terminándose el juego de la vida. Por ejemplo: Figura 1: John Conway; Fuente: lifepatterns. 6 En el caso A B Y C, al no existir una vecindad frontal, el programa no interpreta una conexión dando por tanto el juego terminado. Sin embargo, en los casos D y mi, si existe vecindad próxima, por lo que pueden darse dos casos: muerte o crecimiento, dependiendo en cualquier caso de algoritmo en este caso. Debemos tener en cuenta este tipo de juego como un ejemplo didáctico del funcionamiento y simulación de un Autómata Celular, y no todos los movimientos del ente se ajustan específicamente a los fenómenos del mundo real. 2.1.2 - La importancia de los fractales. Antes de todo, debemos preguntarnos desde un primer momento: ¿Qué es un fractal?, ¿Qué tiene que ver con la Geografía o con los Autómatas Celulares? Los fractales son formas irregulares donde su geometría se ve versada en función de la escala de visión. La irregularidad de sus formas cambia según la distancia desde la cual se observe. El término fue acuñado por Benoit Mandelbrot en 1975. En la actualidad existen numerosos tipos y formas de fractales de los cuales el más popular es el conjunto de Mandelbrot, también denominado: El triángulo Sierpinski (Fig. 2 y 3). Figuras 2 y 3: John Conway; Fuente: lifepatterns. 7 También podemos destacar la alfombra de Sierpinski (Fig. 4): Figura 4: Alfombra de Sierpinski; Fuente: Lifepatterns. Se tratan de formas muy sencillas de construir. En el primer caso, nos basamos en dibujar un triángulo. A su vez dentro de este se vuelve a dibujar sobre su espacio, otro triangulo y así sucesivamente, llegando al infinito computacionalmente. Conforme aumentamos la escala, observaríamos como los triángulos se repiten una y otra vez sin finalización posible. Una parte importante en los fractales es su autosimilaridad, consistente en la propiedad por la que el objeto se recrea así mismo como ocurre con las formas de Mandelbrot. Dentro de la Geografía, los fractales son un modo de uso caracterizado para áreas irregulares, como los litorales o las estructuras urbanas. En el caso de las líneas costeras, conforme se aumenta la escala, las líneas litorales cambian al igual que los fractales, por lo que es un método de uso principal para las costas escarpadas. La línea litoral se diferencia en función de la escala entre lo irregular a escalas más finas y regulares a escalas más distantes. A este hecho costero hay un gran consenso al cual se ha denominado la paradoja de Steinhaus (Nysteun 1966), a pesar de que el asunto fue iniciado por Richardson en 1961, quien describió este fenómeno. En nuestro caso, la aplicación geográfica es la basada en la cuestión urbana. En lo correspondiente a la idea del litoral, ocurre un efecto similar en los límites urbanos de las ciudades. Perkal (1958), realizó mediciones sobre la ciudad de Wroclaw, siendo el uso de 8 los fractales para las zonas de frontera. A pesar de que existen muchos moldeadores de la ciudad y su planificación, el uso de los fractales es útil para los usos del suelo a la hora de conocer las áreas de curvatura similar en los bordes y para comparar ciudades entre sí. 3 – METODOLOGIA 3.1 Métodos de aplicación del Autómata Celular. El uso de entes celulares como simulación del crecimiento urbano debe ser claramente enfocado sobre una realidad como es, en nuestro caso, Madrid. Por ello, para realizar un correcto funcionamiento de los Autómatas Celulares aplicado sobre nuestra área escogida, debemos tener en cuenta como se ha comentado anteriormente, una serie de parámetros o factores como son: 3.1.1 – Accesibilidad. Se utiliza como factor a tener en cuenta sobre la red de transporte (Fig.5). Es de vital importancia sobre todo para conocer si los diferentes usos poseen un potencial de accesibilidad óptimo a la red viaria como, por ejemplo, los usos industriales o residenciales. Para generar un mapa de accesibilidad, se utilizará la aplicación de la Distancia Euclidiana a través de un software GIS. En un mapa de este tipo la distancia es importarte en cuanto a su evaluación. Es por ello, por lo que la cercanía a carreteras tendrá un alto valor. La ecuación matemática para calcular la accesibilidad a partir de la distancia Euclidiana es la siguiente (White et al. 1997): D es la Distancia Euclidiana de cada píxel al más cercano, mientras δj es el coeficiente que nos indica la importancia de la accesibilidad. Por tanto, este valor controla si su evaluación sea mayor o menos, dando una atracción mayor o menos a la red de carreteras en cuanto a su accesibilidad. El conjunto de estos factores tiene la misión de calibrar el modelo sobre el que realizaremos los cálculos, por lo que es de vital importancia su práctica con el fin de conseguir una mayor realidad urbanístico. 9 Figura 5: Mapa de carreteras utilizadas como factor de accesibilidad; Fuente: propia 10 3.1.2 - Factor aleatorio. Se basa en la posibilidad de que una ciudad cuando se desarrolla posee un cierto desorden y orden al mismo tiempo como consecuencia de las diferentes políticas de ordenación desarrolladas durante décadas. Por tanto, el factor estocástico es fundamental en cualquier caso urbanístico (Fig.6). El factor aleatorio (v), aplicado en el estudio es proporcionado por una función de distribución de Weibull como se muestra en la siguiente ecuación (White & Engelen 1993). v = 1+(-ln(1-random)) ·exp(α) 𝛼 es un parámetro escalable que controla el efecto estocástico. Random, se refiere a un número pseudoaleatorio de la distribución uniforme. 3.1.2.1 - Establecer los diferentes tipos de uso del suelo sobre el área de estudio: Con este fin podremos acometer qué atributos queremos analizar a través de su crecimiento. Debemos conocer que elementos son los adecuados para el estudio de un área concreta. 3.1.2.2 - La escala: Este factor es relativamente importante a la hora de realizar los cálculos de estudio, puesto que, si la celda de resolución es alta, las correcciones no van a ser óptimas o con mayor error de distancia debido a la resolución. No es lo mismo en un área de 20 Km2 una red de 50x50, que una de 3x3. La realidad a la hora de escoger una escala es tener en cuenta no solo las capacidades del sistema, sino también de los cálculos o la computación del ordenador que se utilice por lo que, a nivel de eficiencia, los más lógico es utilizar una escala intermedia que nos permita un mayor nivel de cálculo por segundo, así como una malla de resolución acorde al área. En nuestro caso se usará una resolución de 100 x 100. 11 Figura 6: Factor aleatorio sobre la CAM aplicando la herramienta Random Raster para la generación del efecto estocástico; Fuente: Propia. 12 3.1.3 – Vecindad. Es un factor muy importante para los diferentes usos del suelo por diferentes motivos. El más importante es que, si tenemos un uso de suelo de tipo industrial, su crecimiento es mayor o con más atracción sobre aquellos puntos similares, por lo que esas áreas podrán tener una mayor especialización de ese uso durante las generaciones que se active el Autómata Celular. Así mismo, la vecindad se aplica sobre un área, donde el centro posee una mayor atracción que las áreas limítrofes, pues el factor de la distancia limita su atracción a efectos. La vecindad se calcula a partir de la siguiente ecuación: 𝑁𝑖𝑗 𝑡+1 = ∑∑𝑤𝑘𝑥𝑑 𝑡 Donde 𝑤kxd es el parámetro de peso aplicado en un mapa de usos del suelo, t al uso del suelo, k en una posición x, en una zona de distancia d de la vecindad y 𝑁𝑖𝑗 𝑡+1, el factor de transición potencial de la vecindad. 3.1.3.1 – Representación de los usos urbanos en el modelo: Se utiliza como un valor entre 0 y 1 (valor binario) donde se le permite a la celda tener un determinado tipo o uso del suelo en función del área que se analice. En nuestro caso vamos a diferenciar entre suelo urbano (1) y no urbano (0). Todos los valores que no sean 1 y 0 se clasifican como No Data (Fig.7). 13 Figura 7: Cambios en los suelos urbanos en el periodo del año 2000 al 2009; Fuente: Propia. 14 Madrid Land Use es un mapa detallado de los usos del suelo de la CAM, con datos basados en la foto- interpretación con una gran base de datos geográfica sobre el conjunto de datos del área urbana de los años 2000-2006-2009. Su creación se produjo como consecuencia de los grandes cambios en el suelo de la comunidad en el primer período del siglo XXI (Díaz-Pacheco & García-Palomares (2014)). A partir de una capa de usos de la CAM como es Madrid Land Use (MLU), extraemos los valores de suelo urbano y no urbanos mediante una operación de reclasificación. Puesto que en MLU solo están cartografiados los suelos urbanos, el trabajo de reclasificación es sencilla – asignamos el valor 1 a todos los polígonos que tienen un uso asignado, y el valor de 0 a todos los polígonos designados “sin clasificar”. A las zonas verdes urbanas, como la Casa de Campo, asignamos un valor de “NA”, eliminándolas de la modelización. De este modo, suponemos que todo puede cambiar, menos las zonas verdes urbanas, y que todo suelo que no es ya urbano, está libre para ser urbanizado. En algunas ciudades, esto no sería una suposición adecuada. Sin embargo, Madrid se encuentra rodeado de terrenos agrícolas , donde (en la mayoría de los casos) se ha permitido construir sin demasiadas restricciones. 3.1.4 – Aptitud. El último factor que se incluye es la aptitud física del terreno que actúa como factor limitante sobre el desarrollo urbano, afectando a la posibilidad de que un uso del suelo pueda establecerse en un determinado lugar. En nuestro caso, a partir de un Modelo Digital de Elevaciones (DEM), se extraerá un mapa de pendientes para conocer el ángulo de esta (Fig.8 y9 ). La pendiente será reclasificada según un valor del 1 al 10 siendo 1 alta pendiente y 10 poca pendiente(Clarke et al. 1997; Schneider et al. 2001). En cuanto a la reclasificación de la pendiente, se hace teniendo en cuenta el parámetro por el cual, el mayor valor se le asigna a la zona donde mejor posición de construcción de suelo urbano se considere. En este caso la sierra, con una pendiente de valor 1 sería un lugar prácticamente difícil, al no ser un suelo llano. 15 Figuras 8 y 9: Mapa de altitud (en metros) y pendientes de la comunidad de Madrid en 2D y 3D; Fuente: propia. 16 Mapa de Pendiente Mapa de Altitud Valores de la pendiente. 9: Poca pendiente. 1: Mucha pendiente. Valores de altitud. Expresadas en metros. Modelado 3D de la Comunidad de Madrid. 17 3.1.5 - Cálculo del potencial de transición. Mediante la accesibilidad (A), vecindad (N), aptitud (S), zonificación (Z, es opcional en algunos casos) y aleatoriedad (v) calculados en anteriores pasos, se determina este factor a partir de la siguiente ecuación: 𝑃𝑖𝑗 𝑡+1 = 𝑓(𝑁𝑖𝑗 𝑡 , 𝐴𝑖𝑗 𝑡 , 𝑆𝑖𝑗 𝑡 , 𝑣𝑖𝑗 𝑡 ) La función (f) que se aplica dependerá de la importancia de los factores en el caso de estudio concreto, En este trabajo, se experimentó con dos funciones de transición distintas: 1) En la fase de calibración, se multiplicó cada factor o parámetro sin normalizar a través de esta secuencia en R: model_TTP <-(model_accessibility+1) * (model_suitability+1) * (model_nhood+1) * (model_random) 2) Para las simulaciones finales, se aplicó la suma de los factores normalizados: (Factor de vecindad) topn <- model_nhood-(minValue(model_nhood)) bottomn <- maxValue(model_nhood)-minValue(model_nhood) TPn <- (topn/bottomn) (Factor aleatorio) topr <- model_random-(minValue(model_random)) bottomr <- maxValue(model_random)-minValue(model_random) TPr <- (topr/bottomr) (Factor de accesibilidad) topa <- model_accessibility-(minValue(model_accessibility)) bottoma <- maxValue(model_accessibility)-minValue(model_accessibility) TPa <- (topa/bottoma) (Factor de aptitud) 18 tops <- model_suitability-(minValue(model_suitability)) bottoms <- maxValue(model_suitability)-minValue(model_suitability) TPs <- (tops/bottoms) (Suma de los diferentes factores) model_TTP <- CS + TPn + TPr + TPa + TPs La finalidad de haber realizado este cambio es obtener en los mapas estáticos y dinámicos finales, una serie de cambios visibles, para reconocer en la simulación las diferencias entre ambos modelos. La vecindad en muchos casos, como se ha podido comprobar en los mapas de calibración (en el apartado 4.2.1.), contiene una mayor influencia en numerosas áreas del mapa de la CAM, provocando que la accesibilidad no fuera los suficientemente fuerte como para tener una presencia óptima sobre la versión final. Es por ello por lo cual se ha decidido variar el factor de transición potencial únicamente en los modelos finales al ser los mapas con la calibración previamente calculada y establecida. Siguiendo esto, se podrá conocer la evolución del estado de cada píxel o celda determinado, gracias a este factor. Como último paso, se debe efectuar un bucle al sistema con el fin de que los parámetros establecidos entre diferentes fechas surjan efecto, es decir, se observe la evolución en el tiempo del espacio urbano cambiante sobre la cartografía, en este caso, de la Comunidad de Madrid. 4 - SIMLANDER SIMLANDER es un software desarrollado principalmente como simulación de usos de suelo, con la finalidad de conocer, mediante software R, los cambios de los usos del suelo mediante una serie de cálculos y calibraciones realizadas por el usuario (Hewitt. et al. 2013). En nuestro caso, se calibrará el sistema con el fin de simular el suelo urbano a desarrollar entre los años 2000 a 2009 con una accesibilidad, primero estática y otra dinámica. Una vez se tenga el sistema calibrado, se realizará una simulación mayor desde el año 2000 al 2050 para conocer de esta forma, un desarrollo del suelo urbano lo más real posible. 19 4.1 - El uso de Software R como GIS. El uso de un lenguaje de programación como es R aplicado a los Sistemas de Información Geográfica nace con la premisa de establecer nuevos métodos de cartografía. R como software aplicado a los SIG es usado para la geoestadística, análisis espacial, tratamiento de datos en los formatos que se presente , modificarlos, etc. En nuestro caso de aplicación de este software, se basa en la idea de facilitar la manera de tratar los datos y modificar los cálculos de una manera más sencilla que con software ya desarrollado mediante Apis o scripts, sin la posibilidad de modificarlos a nuestro gusto. Con todo esto y con las herramientas necesarias en R, se ha creado un script totalmente acomodado a nuestra tarea (SIMLANDER 4.2.1), en la cual se pueden cambiar los parámetros de una manera más ágil sin la necesidad de repetir los procesos una y otra vez en cada caso que se necesite. Con la creación de variables y parámetros, la tarea se realiza más sencilla al modificar en el código simplemente las variables asignadas, lo que se traduce en mayor optimización y celeridad. Por tanto, el uso de R como un Sistema de Información Geográfica, no es más que otra forma de software, en este caso Open Source, que nos permite realizar las mismas tareas que un programa creado como QGIS o ArcGIS, pero, con la diferencia de que en R nosotros somos los que decidimos como hacer todo sin la necesidad de seguir unas pautas específicas. Nosotros somos los encargados de cómo obtener los resultados finales. 4.2-Aplicación de SIMLANDER mediante R (calibración del sistema). En la calibración de SIMLANDER para la simulación del suelo de la Comunidad de Madrid, es necesario seguir una serie de pasos previos, comenzando por la limpieza de datos que se usarán evitando en la medida de los posible , tener solo lo justo y necesario para evitar errores de cálculo. Tras la aplicación de los diferentes procesos explicados anteriormente en el apartado de metodología, se debe comenzar con la misión de calibrar la accesibilidad y la vecindad. En el caso de vecindad, se comenzará usando un programa en Excel (creado por Richard Hewitt, Fig.10) para conocer, primeramente, que tipo de vecindad se quiere utilizar, probando diferentes cálculos, basándose en tener en cuenta la influencia en función de la distancia. A partir de la gráfica se puede observar, como la influencia va disminuyendo en función de la distancia que nosotros mismos propongamos, aunque en este caso la utilizada, se usa una influencia mayor desde los 0 metros a los 100 metros. Sin embargo, a partir de esta distancia y en aumento de esta, la influencia baja drásticamente hasta 0 a los 300 metros de distancia. En nuestro caso se aplicarán 3 vecindades diferentes que se comentarán en el apartado 4.2.1.1 en conjunción con la accesibilidad y los mapas calibrados. 20 Figura 10 :Aplicación de Excel para el cálculo de la vecindad, pudiendo establecer el área de influencia en metros; Fuente propia. En el caso de la accesibilidad, se han utilizado 3 cálculos diferentes en el Script: Accesibility_1 <- (1+(distroad/1))^-1 Accesibility_2 <- (1+(distroad/1000))^-1 Accesibility_3 <- (1+(distroad/1000000))^-1 La diferencia principal entre los diferentes métodos aplicados sobre la accesibilidad, radican en la forma de la curva de la caída con la distancia, siendo o bien muy brusco (accesibilidad 1, atracción muy fuerte al lado de las carreteras, muy débil a poca distancia de ellas), o bien mucho más suave (accesibilidad 3, atracción ligeramente más fuerte al lado de las carreteras, pero aún importante a gran distancia de ellas). Se experimentó con 21 estos tres tipos de accesibilidad, con el fin de conocer el comportamiento de Autómata Celular en diferentes casos. Por tanto, a continuación, se generarán a través del script de SIMLANDER, una serie de mapas con los diferentes tipos de vecindad y accesibilidad. En el apartado ANEXOS, mostramos el script de SIMLANDER usado en R. 4.2.1 Accesibilidad estática (2000-2009). Tras obtener los diferentes parámetros a los que alimentar al script con sus respectivos factores, acudiremos a R a correrlo con las diferentes variaciones de accesibilidad y vecindad para conocer cómo se desarrollará el mapa de la CAM entre los años 2000-2009. Así mediante el uso de estos parámetros escogidos, se procederá a realizar sucesivas calibraciones del modelo para obtener el más acorde a la simulación real del año 2009. Con las calibraciones obtenidas a continuación, procederemos a utilizar Map Comparison Kit (Visser and de Nijs 2006), un software gratuito, disponible en: http://mck.riks.nl/, para conocer las diferencias entre los distintos mapas generados, De esta manera, una vez se presente la mejor calibración, se realizará a posterior un estudio a futuro sobre el desarrollo del suelo urbano hasta el año 2050. Las comparaciones de hicieron de dos maneras distintas, mediante inspección visual, ayudado con mapas de colores, y mediante un indice de similitud llamado “Kappa Simulation” (Ksim) (Van Vliet et al 2011). La figura 11, nos mostrará el uso de la accesibilidad 1 con las vecindades 1, 2 y 3. La figura 12, nos mostrará el uso de la accesibilidad 2 con las vecindades 1, 2 y 3. La figura 13, nos mostrará el uso de la accesibilidad 3 con las vecindades 1, 2 y 3. La figura 14, nos mostrará el uso de la accesibilidad 4 con las vecindades 1, 2 y 3. http://mck.riks.nl/ 22 Accesibilidad 1 → W1 → W2 → W3 Figura 11: mapas de calibración de accesibilidad 1 y vecindades 1,2 y 3; Fuente: propia. 23 Accesibilidad 2 → W1 → W2 → W3 Figura 12: mapas de calibración de accesibilidad 2 y vecindades 1,2 y 3; Fuente: propia. 24 Accesibilidad 3 → W1 → W2 → W3 Figura 13: mapas de calibración de accesibilidad 3 y vecindades 1,2 y 3; Fuente: propia. 25 Accesibilidad 4 → W1 → W2 → W3 Figura 14: mapas de calibración de accesibilidad 4 y vecindades 1,2 y 3; Fuente: propia. 26 4.2.1.1 - Comparativa de los mapas estáticos entre los años 2000-2009. Resultados de la aplicación de la estadística de Kappa Simulation (Ksim), de Map Comparison Kit, a través de las celdas raster de cada uno de los mapas creados. Cuanto más alto sea este valor, mayor coincidencia hay entre el mapa generado a través de la simulación del año 2009 y el mapa real de referencia para el año 2009. Tabla 1: Resultado de los mapas a través de Map Comparison Kit; fuente propia. Para comprender los mapas es necesario conocer que: -Las celdas en color rojo representan áreas de Madrid en el año 2009 real. -Las celdas en color azul representan áreas de Madrid de 2009 simuladas y creadas por el Autómata Celular. -Las celdas en color verde representan las áreas de Madrid que se encuentran en ambos casos anteriores (rojo y azul), tanto la real como la simulada. Podemos afirmar que el mapa más aproximado al crecimiento real de la CAM es el correspondiente con la Vecindad 2 y un valor de Accesibilidad 1 (Fig.11, mapa 2 ). Sin embargo, se pueden extraer diferentes parámetros a comentar sobre el conjunto de mapas en función del cálculo de la vecindad o accesibilidad usadas en cada caso para conocer que zonas o lugares se pueden haber visto más afectadas por crecimiento del suelo o, en caso contrario, por la reducción de este. En el caso de los mapas con Vecindad 1 siendo su cálculo aplicado: Valor de Ksim Vecindad 1 Vecindad 2 Vecindad 3 Accesibilidad 1 0.010 0.136 0.134 Accesibilidad 2 0.012 0.105 0.115 Accesibilidad 3 0.010 0.135 0.135 Accesibilidad 4 0.013 0.075 0.084 27 (0.000,0.000,0.005,0.000,0.000, 0.000,0.031,0.050,0.031,0.000, 0.005,0.050,0.500,0.050,0.005, 0.000,0.031,0.050,0.031,0.000, 0.000,0.000,0.005,0.000,0.000), nr=5,nc=5) Observamos que en función de la accesibilidad que apliquemos, vemos una mayor o menor concentración. Por tanto, cuando se aplica la Vecindad 1 sobre las diferentes accesibilidades usadas, la dispersión es alta, en cualquier caso. Esto se debe a que la proporción de la distancia a los principales viales se hace de manera más desproporcionada, debido a que la vecindad 1, sobre la matriz creada en función de la ley de Tobler, no se le ha dado una fuerza lo suficientemente alta para que la atracción de las nuevas celdas de uso urbano en los núcleos sea suficiente para compensar la atracción hacia las vías de transporte. A su vez, las vías de transporte por su accesibilidad actúan en función de la matriz de vecindad 1 (en este caso), que como hemos dicho y diferenciándola de las otras vecindades, le otorga un valor de Ksim realmente bajo y nada realista a las calibraciones posteriores que veremos. Con el uso de la Vecindad 2, los cálculos y la simulación son más acordes a la realidad de crecimiento sobre la CAM. En este caso los valores de vecindad 2 utilizados han sido: (0,000,0.000,0.500,0.000,0.000, 0.000,3.136,5.000,3.136,0.000, 0.500,5.000,50.000,5.000,0.500, 0.000,3.136,5.000,3.136,0.000, 0.000,0.000,0.500,0.000,0.000), nr=5,nc=5) A diferencia de la vecindad 1, la influencia entre las celdas es mayor entre ellas, sin embargo, no lo suficientemente diferenciada. Por tanto, la influencia entre las celdas adyacentes es mayor al caso 1. 28 Usando la accesibilidad 2, el patrón que se genera de la simulación resulta menos realista al crecimiento en los 9 años generados, y existen diferencias importantes con respecto a los mapas de accesibilidad 1: los centros de la ciudad de Madrid. Podemos ver como en el caso de la vecindad 1, el centro de Madrid se encuentra en color verde con la vecindad 2 y 3, mientras que en los mapas de accesibilidad 2 (Fig.10), con vecindad 2 y 3 (Fig. 16), los centros no se encuentran fielmente representados. Figura 15. Accesibilidad 1, vecindad 2 y 3; Fuente: propia. Figura 16: Accesibilidad 2, vecindad 2 y 3; Fuente: propia. Efectivamente observando los mapas, se ve claramente como el centro de la ciudad de Madrid se encuentra muy poco representado en la figura 16. Esto se debe por el hecho de que el aumento de la accesibilidad a las carreteras radiales propicia una mayor influencia 29 sobre los viales que sobre el propio centro debido a que estos viales, no entran en el centro de Madrid desde el Km 0, sino sobre los entre 5- 10 Km de radio tomando como referencia el centro de la Puerta del Sol. Por tanto, a efectos mayores lo que ocurre, es que la influencia se desplaza hacia los viales de accesibilidad debido que en los mapas de accesibilidad 2 (Fig.12), la fórmula aplicada sobre la accesibilidad es 100 veces mayor a la usada en los mapas de accesibilidad y, por ende, la influencia tiende a aumentar sobre las radiales primordialmente. Si ahora nos vamos a ver los resultados con la Vecindad 3, se explicará la matriz utilizada en este caso siendo: (0.000,0.000,0.500,0.000,0.000, 0.000,29.496,50.000,29.496,0.000, 0.500,50.000,500.000,50.000,0.500, 0.000,29.496,50.000,29.496,0.000, 0.000,0.000,0.500,0.000,0.000), nr=5,nc=5) La matriz se encuentra con una misma área de influencia que las anteriores, con la diferencia de que la fuerza entre las celdas es mayor. A pesar de que las fuerzas de la celda sean mayores, los aledaños no tienen una mayor influencia que las anteriores debido a que esta misma solo llega hasta los 300 m máximo, por lo que los mapas van a ser similares a los anteriores, posiblemente muy comparables con los de vecindad 2. Sin embargo, la vecindad 3 en los mapas creados, dan un nivel de Ksim superior a los de accesibilidad 1 (Fig.11) y algo más próximos a los de vecindad 2. La similitud con la vecindad 2 se basa en que los números o datos empleados sobre la matriz son muy similares, simplemente añadiéndose mayor o menor fuerza de atracción sobre la misma área de celdas, por lo que no se ve incrementada en cuanto a distancia. Si se modificara la distancia de atracción, si se producirá un cambio en el mapa ampliándose el radio. Como caso final, también se ha calculado una accesibilidad 4 (Fig.14) caracterizada por ser su fórmula aplicada a restringir ciertos datos: con( (distroad < 300), 1, con( (300 <= distroad) & (distroad < 1500), (0.5*(1 + cos((3.14159*(distroad-300))/1200))),0)) 30 A través de esta fórmula aplicada, lo que se realiza son condicionales de distancia, con un intervalo de 300 m como se ha aplicado en la vecindad en Excel (Fig.10). Por tanto, con este tipo de accesibilidad que se ha realizado en los mapas de accesibilidad 4 (Fig.14), se representa algo más real la simulación del Autómata Celular a la vez que le estamos comunicando una serie de condiciones sobre donde debe actuar a tener en cuenta con la influencia de la vecindad y la accesibilidad misma. A pesar de que los datos reales de Ksim con esta fórmula no son lo más altos, si se debe tener en cuenta que la simulación es muy realista, por lo que una calibración sobre la accesibilidad 4 (Fig.14) y aumentáramos el radio de influencia de la vecindad, así como un aumento generalizado de los viales, el mapa crecería. Pero se debe tener en cuenta, que, en este caso, solo tomamos como viales importantes las radiales y las carreteras de circunvalación M-30, M-40, M-45 y M-50. 4.3 - Modelización dinámica y estática entre los años 2009-2050. Tras las calibraciones en el apartado anterior, se realizará a continuación dos mapas finales basándonos en la posibilidad de modificar la accesibilidad a través de un modelo dinámico en el cual, por cada año que el bucle del script realice la transformación cartográfica, tenga en cuanto el cambio de accesibilidad generada. Por tanto, se comparará el modelo estático utilizado en la calibración del modelo y otro modelo dinámico, para los años 2009-2050 conociendo la posibilidad de crecimiento del suelo urbano sobre la CAM. En ambos casos se mantuvieron las mismas reglas de vecindad, obtenidas con la calibración que mejor puntuación daba según la estadística de Ksim. ¿Qué diferencia el modelo estático al dinámico? La principal diferencia entre ambos modelos se basa en los cambios de accesibilidad progresivos que se suceden en el modelo dinámico, algo que no ocurre con el estático. Es decir, si por un lado vemos como se desarrolla la accesibilidad en el modelo estático, ateniéndonos al script: model_accessibility <- accurb acc<-calc(accessibility, funacc) end calculate accessibility Para el modelo estático, el factor de accesibilidad solo lo realiza como uno de los pasos previos al inicio del script, para el desarrollo de los factores que luego se usaran para el factor de transición: 31 model_TTP <- (model_accessibility+1)*(model_suitability+1)*(model_nhood+1)*(model_random) Como hemos comentado anteriormente, al tener una accesibilidad estática sobre las vías de circulación, la accesibilidad vas a ser exactamente la misma en cada uno de los procesos por los cuales el Autómata Celular se desarrolla con los progresivos años dados (2009-2050). Por otro lado, si vemos en el script el modelo dinámico y observamos como se realiza la accesibilidad, directamente acudimos al bucle donde se van generando los mapas por año: #begin dynamic accessibility block print("calculating dynamic accessibility.......") model_accessibility <- accurb/i d <- gridDistance(Tn, origin = 1 d2 <- abs(d-(cellStats(d,max))) d3 <- d2/(cellStats(d2,max)) d4 <- model_accessibility + (d3/1000) d5 <- d4/(cellStats(d4,max)) model_accessibility <- d5 A diferencia del modelo estático, la accesibilidad, se calcula en el momento en el que se genera el mapa año por año, por lo que, en este caso, se tiene en cuenta los cambios de vecindad y accesibilidad a medida que se corre el modelo durante los sucesivos años. Teniendo en cuenta esto, la simulación tiende a ser más real usando el modelo dinámico por esto mismo; la posibilidad de que a medida que se desarrolla el Autómata año a año, el suelo urbano se va modificando a razón de las reglas establecidas, lo que conlleva a que la vecindad al calcularse por cada cambio en el suelo sea más real y acorde a la vecindad establecida con el crecimiento. Aun con todo, se observarán a continuación que cambios existen entre ambos modelos a modo de comparativa a través de las aplicaciones de Map Comparison Kit y Fragstats (McGarigal 1995). 32 Modelo Estático Figura 17: Mapa de la CAM (2009-2050) con el modelo estático aplicado en R; Fuente: propia. 33 Modelo Dinámico Figura 18: Mapa de la CAM (2009-2050) con el modelo dinámico aplicado en R; Fuente: propia. 34 Figura 19: Mapa final sobre la simulación con AC sobre la comunidad de Madrid. Año 2050; Fuente: propia. 35 Figura 20: Mapa final sobre la simulación con AC sobre la comunidad de Madrid. Año 2050; Fuente: propia. 36 4.3.1 - Comparación con Map Comparison Kit. ¿Como podemos saber si un modelo es mejor o se asemeja más a la realidad de la CAM? En este caso, es complicado dictaminar si las reglas del Autómata Celular son perfectas, pues, a pesar de que el modelo escogido en la calibración sea bueno, siempre pueden crearse mejoras en el script teniendo en cuenta numerosas variables más que generaran unos cambios más visibles. En nuestro caso si acudimos a observar las diferencias entre ambos modelos a través de MCK, existen cambios en los usos del suelo entre ambos modelos. Antes de todo se debe conocer qué tipo de mapa es de cada color para conocer las diferencias: En color verde, son las áreas igualitarias en ambos modelos, tanto estático como dinámico. En color rojo, son las celdas del modelo dinámico. En color azul, las celdas del modelo estático. Echando un primer vistazo, se observan cambios más puntuales en toda la franja comprendida entre el nordeste hasta el sur de la comunidad, pero no son los únicos casos. Esto responde al hecho de que, en estas áreas, al tener también un tipo de relieve más llano que la zona Norte y Oeste, propicia un mayor desarrollo de suelo urbano, así como también la posibilidad de que pueblos que ahora en el 2020, comienzan a aumentar de tamaño, se consideren ciudades, o incluso como se puede ver en la zona del corredor del Henares, la unión completa entre los municipios de Torrejón de Ardoz, Alcalá de Henares, San Fernando de henares, Coslada. Figura 21: Corredor del Henares con la superposición de los modelos; Fuente: propia. 37 En la zona de las Vegas, Arganda del Rey, Nuevo Baztán, etc. situado al Este y el Sureste, existe una gran concentración de pequeños núcleos de población que en la actualidad se encuentran en desarrollo por la creación de ensanches con viviendas, en la mayoría de los casos unifamiliares al estilo americano. Es una zona interesante de observar y posiblemente la más expuesta al desarrollo. Figura 22: Área este de la CAM, en superposición con los modelos. Fuente: propia. Sobre la zona sur de la CAM, el mayor cambio se observa en los municipios del extremo sur como son Parla, Pinto, Valdemoro o incluso Torrejón de la Calzada. Figura 23: Área sur de la cam, en superposición con los modelos; Fuente: propia. 38 Otro de los núcleos interesantes de crecimiento que, si ocurre sobre la zona Norte, es la correspondiente con los municipios de Norte a Sur desde San Sebastián de los Reyes hasta Talamanca del Jarama. En este caso la accesibilidad y la posibilidad de existencia de suelo disponible para ser urbanizado en mayor a pesar de que en esta área comienza un aumento progresivo de altitud e irregularidad del terreno a partir de Talamanca. Figura 24: Área norte de la CAM, orientada en horizontal. De izquierda a derecha, San Sebastián de los Reyes- Talamanca del Jarama; Fuente: propia. Y finalmente, no menos importante, aunque con menos evolución de suelo, el área oeste y norte, los núcleos que más cambian son los correspondientes con Collado Villalba, al Norte, la zona de Soto del Real y Manzanares y también destacar los municipios situados al suroeste de la CAM como es el entorno de Villa del Prado. Figura 25: Área noroeste de la CAM, con los núcleos de Collado Villalba y Soto del Real al norte. Al sur Villa del Prado; Fuente: propia. 39 Figura 26: Mapa de la CAM contrapuesto en Map Comparison Kit con los modelos dinámico y estático; Fuente: propia. Mapa de la comunidad de Madrid a través de MCK 40 Como conclusiones sobre la comparativa, tanto el modelo dinámico como estático son muy similares, en cuanto al crecimiento urbano. A pesar de que existan diferencias entre ambos modelos, se puede observar cómo los dos crecen de forma muy similar en los lugares que se han observado. Por tanto, en estos casos expuestos, la opción de escoger entre un mapa estático (Fig.17 y 19) o uno dinámico (Fig.18 y 20), conllevaría a realizar modificaciones sobre ambos modelos, pues de esta manera el aumento del suelo urbano en ambos casos es muy similar. Esta similitud se atiende al peso que existe en la vecindad en los usos del suelo (igualitarios en ambos modelos). De hecho, si se fija un mapa de los viales en formato vectorial, se vería como la accesibilidad en el modelo dinámico, posee una mayor atracción, por sus cambios. Por el contrario, en el estático la accesibilidad permanece invariable. Si se tuviera que escoger un modelo, siempre va a ser más propicio el uso de una accesibilidad dinámica por el hecho de que este tipo de modelo se retroalimenta del mapa que va generando en cada sucesión a partir de los cálculos de vecindad y accesibilidad en cada año, por lo que se aplica una nueva configuración en el suelo urbano. En el caso del modelo estático solo actúa como factor dinamizante, la vecindad. Esta es la principal diferencia entre ambos modelos. 4.3.2 - Comparación con Fragstats. Fragstats permite conocer si el resultado obtenido con los modelos estático y dinámico, poseen ciertas similitudes sobre el patrón espacial con respecto a sí mismo al mapa de 2009 como comparación. (MCGARIGAL, K. (1995)). El motivo de probar este tipo de comparación y no utilizar el índice de Ksim como en la fase de calibración, se debe a que Ksim requiere un mapa de referencia acorde con la fecha simulada, que era 2009 en nuestro caso. Lógicamente, con una simulación hacia el futuro, no contamos con un mapa real de la fecha de referencia (2050), y se tiene que buscar otras técnicas para saber la calidad de nuestra simulación. En caso de ser similares, el resultado sería óptimo, demostrando que la innovación con la accesibilidad dinámica no cambia el patrón ni la estructura de la simulación, que, a su vez, es muy similar al patrón real observado en lu2. A través de Fragstats, se compararon los 3 mapas en formato GeoTIFF, para conocer el valor de cada uno de los mapas siendo los valores: Tabla 2 con los resultados del modelo real 2009 en fragstats; Fuente: propia. 41 Tabla 3 con los resultados del modelo estático (2050) en Fragstats; Fuente: propia. Tabla 4 con los resultados del modelo dinámico (2050) en Fragstats; Fuente: propia. Se observa un patrón fractal muy similar entre los mapas de las tablas 2, 3 y 4. Existen ciertas diferencias, sin embargo, tal y como podemos observar, el valor de la dimensión fractal es de 1.02 en todos los casos. Tal y como nos indican los autores: La investigación sobre la naturaleza fractal de las ciudades es anterior a los modelos de crecimiento urbano de CA y varios autores recientes la han empleado para medir el grado de ajuste en el modelado de simulación del uso del suelo (Hewitt y Díaz-Pacheco, 2017; Newland et al., 2015). La dimensión fractal masiva mide el grado de linealidad de la clase urbana en general, en la que los objetos que llenan el plano como círculos o cuadrados tienen un valor de 2.0 y una línea tiene un valor de 1.0. Por tanto, nuestro caso de estudio entre ambos modelos comparados con la simulación real de 2009, si poseen bastante similitudes y por ende, son modelos completamente válidos a nivel fractal y estadístico. 42 5 - CONCLUSIONES Los Autómatas Celulares aplicados como método de análisis geográfico, proyecta grandes posibilidades que se han comprobado a través de este estudio. Se ha analizado como el suelo urbano de la Comunidad de Madrid variará hasta 2050, a través de dos modelos distintos. En el primer modelo, aunque las normas de vecindad se aplican en cada paso de la simulación, basándose en la configuración de los usos del suelo del paso anterior, el factor de accesibilidad se mantiene estático, es decir, sin cambiar entre el primer paso temporal de la simulación y el ultimo. El segundo modelo, tanto las reglas de vecindad como el factor de accesibilidad son dinámicas, ya que este último se ha introducido dentro del bucle temporal, haciendo que la accesibilidad se actualice en cada paso de la simulación. El objetivo de haber realizado estos dos modelos es comprobar que efectivamente, ambos modelos son completamente válidos a la hora de realizar este tipo de estudios basados en el crecimiento urbano. A través de un análisis del patrón espacial de las simulaciones mediante el software Fragstats, se ha determinado que el valor de los mapas en cuanto a su fractalidad, son compatibles y por tanto ambos tienen la posibilidad de adecuarse a los factores de una manera u otra siendo sus resultados válidos. Sin embargo, si se tiene que escoger un modelo más acorde con los cambios, no solo de suelo urbano, (otros tipos de suelo como el vegetal, por ejemplo) el modelo dinámico se considera de mayor realismo al tener en cuenta los cambios de accesibilidad por año y no por fecha en concreto, como ocurre en el estático. Sería necesario generar posiblemente un mayor número de factores o reglas hacia el AC, para hacerlo más realista, por lo que resultaría imprescindible una recalibración. En nuestro caso, la apreciación de ambos modelos no difiere por el hecho de que tienen una serie de factores igualitarios. Si se opta por asignar valores dinámicos sobre otro tipo de usos como el transporte de ferrocarriles o redes de carreteras secundarias (factor de accesibilidad), la simulación probablemente cambiaria drásticamente. Nuestro objeto de estudio ha sido más simple al tener en cuenta la accesibilidad como la red de carreteras principal en primer plano junto al suelo urbano y la vecindad aplicada según una serie de parámetros previamente establecidos. Hemos tenido en cuenta la accesibilidad del suelo urbano hacia los viales, a medida que va creciendo cada año. Sería interesante comprobar a rasgos mayores, la introducción de nuevas variables sobre la Comunidad de Madrid y observar como el Autómata varía por áreas o incluso recrearlo sobre diferentes comarcas de la región y comprobar su comportamiento aislado. El uso del Autómata como parte fundamental del proceso de análisis geográfico resulta muy satisfactorio a la hora de trabajar con el una vez calibrado, puesto que su aplicación, no tiene que ser única y exclusivamente al mundo de la geográfica, sino a otras materias, como es el caso de la Biología en el análisis de movimientos de las especies o su desarrollo sobre un área con el fin de conocer su localización óptima aproximada. 43 6 – REFERENCIAS BIBLIOGRAFICAS 1. B. BHATTA. MODELLING OF URBAN GROWTH BOUNDARY USING GEOINFORMATICS. INTERNATIONAL JOURNAL OF DIGITAL EARTH, 2 (4) (2009), PP. 359-381. 2. C.G.LANGTON (1986). STUDYING ARTIFICIAL LIFE WITH CELULAR AUTOMATA. DEPARTMENT OF COMPUTER AND COMMUNICATION SCIENCE. UNIVERSITY OF MICHIGAN. (PHYSICA 22D PP 120-149). 3. EDWARD N. LORENZ (1962). DETERMINISTIC NON PERIODIC FLOW. MASSACHUSETTS INSTITUTE OF TECHNOLOGY. 4. ENGELEN, G., WHITE, R., ULJEE, I., Y DRAZAN, P. (1995). USING CELLULAR AUTOMATA FOR INTEGRATED MODELLING OF SOCIO-ENVIRONMENTAL SYSTEMS. ENVIRONMENTAL MONITORING AND ASSESSMENT, 34(2), 203-214. 5. HAKEN, HERMANN & PORTUGALI, JUVAL. (2021). CITIES AS HYBRID COMPLEX SYSTEMS. 10.1007/978-3-030-63457-5_2. 6. HEWITT, R., DÍAZ PACHECO, J. AND MOYA GÓMEZ, B., (2013), A CELLULAR AUTOMATA LAND USE MODEL FOR THE R SOFTWARE ENVIRONMENT. 7. JOHN VON NEUMANN & ULAM S. THEORY OF SELF REPRODUCING AUTOMATA (1948). 8. LUGARESARESTI, JOSEBA. (2005). SISTEMAS FRACTALES, CAOS Y HOLÍSTICA EN EL ANÁLISIS TERRITORIAL: LA GEOGRAFÍA BIODINÁMICA. LURRALDE: INVESTIGACIÓN Y ESPACIO, ISSN 0211-5891, Nº 28, 2005, PAGS. 11-30. 9. MAJID SHADMAN ROODPOSHTI, RICHARD J. HEWITT, BRETT A. BRYAN (2020). TOWARDS AUTOMATIC CALIBRATION OF NEIGHBOURHOOD INFLUENCE IN CELLULAR AUTOMATA LAND-USE MODELS. COMPUTERS, ENVIRONMENT AND URBAN SYSTEMS VOLUEM 79, JANUARY 2020, 101416. 10. M BATTY, P.A. LONGLEY.(1986). FRACTAL-BASED DESCRIPTION OF URBAN FORM. ENVIRONMENT AND PLANNING B: PLANNING AND DESIGN, 1987, VOLUME 14, PAGES 123- 134. DEPARTMENT OF TOWN PLANNING, UNIVERSITY OF WALES. INSTITUTE OF SCIENCE AND TECHNOLOGY, CARDIFF, CF1 3EU, WALES. 11. R. WHITE, G. ENGELEN & I ULJEE (1997). THE USE OF CONSTRAINED CELULAR AUTOMATA FOR HIGH-RESOLUTION MODELLING OF URBAN LAND-USE DYNAMICS. ENVIRONMENT AND PLANING. (VOL. 24 PP 323-343). 12. ULAM, S. (1952, APRIL). RANDOM PROCESSES AND TRANSFORMATIONS. PROCEEDINGS OF THE INTERNATIONAL CONGRESS ON MATHEMATICS (VOL. 2, PP. 264-275). PUBLIC SCHOOL PUBLISHING. https://www.sciencedirect.com/science/journal/01989715 44 13. VON NEUMANN, J. (2017). THE GENERAL AND LOGICAL THEORY OF AUTOMATA. IN SYSTEMS RESEARCH FOR BEHAVIORAL SCIENCE SYSTEMS RESEARCH (PP. 97-107). ROUTLEDGE. 14. WHITE, R. AND G. ENGELEN (1993) CELLULAR AUTOMATA AND FRACTAL FORM: A CELLULAR MODELING APPROACH TO THE EVOLUTION OF URBAN LAND-USE PATTERNS, ENVIRONMENT AND PLANNING A, 25: 1175-1199. 15. TOBLER W.R., « CELLULAR GEOGRAPHY », P. 379-386 IN S. GALE, G. OLSSON (DIR.), PHILOSOPHY IN GEOGRAPHY, REIDEL PUBLISHING COMPANY, DORDRECHT, HOLANDA, 1979. 16. MAS, J-F., (2018), ANÁLISIS ESPACIAL CON R: USA R COMO UN SISTEMA DE INFORMACIÓN GEOGRÁFICA, EUROPEAN SCIENTIFIC INSTITUTE, 114 P. 17. DÍAZ-PACHECO, J. & GARCÍA-PALOMARES, J.C. (2014) A HIGHLY DETAILED LAND- USE VECTOR MAP FOR MADRID REGION BASED ON PHOTO-INTERPRETATION, JOURNAL OF MAPS, 10:3, 424-433, DOI: 10.1080/17445647.2014.882798 18. VISSER, H., & DE NIJS, T. (2006). THE MAP COMPARISON KIT. ENVIRONMENTAL MODELLING & SOFTWARE, 21(3), 346-358. 19. VAN VLIET, J., BREGT, A. K., & HAGEN-ZANKER, A. (2011). REVISITING KAPPA TO ACCOUNT FOR CHANGE IN THE ACCURACY ASSESSMENT OF LAND-USE CHANGE MODELS. ECOLOGICAL MODELLING, 222(8), 1367-1375. 20. MCGARIGAL, K. (1995). FRAGSTATS: SPATIAL PATTERN ANALYSIS PROGRAM FOR QUANTIFYING LANDSCAPE STRUCTURE (VOL. 351). US DEPARTMENT OF AGRICULTURE, FOREST SERVICE, PACIFIC NORTHWEST RESEARCH STATION. https://doi.org/10.1080/17445647.2014.882798 45 7-ENLACES WEB HTTPS://SIMLANDER.WORDPRESS.COM/ HTTPS://WWW.UNOCERO.COM/CIENCIA/DE-FRACTALES-Y-AUTOMATAS-CELULARES/ HTTPS://WWW.XATAKACIENCIA.COM/MATEMATICAS/QUE-SON-LOS-FRACTALES-Y-COMO- SE-CONSTRUYEN HTTPS://WWW.IBIBLIO.ORG/LIFEPATTERNS/OCTOBER1970.HTML HTTPS://ES.WIKIPEDIA.ORG/WIKI/JUEGO_DE_LA_VIDA HTTPS://WWW.SCIENCEDIRECT.COM/SCIENCE/ARTICLE/PII/S0198971519300948?VIA%3D HUB HTTPS://WWW.UMASS.EDU/LANDECO/RESEARCH/FRAGSTATS/DOCUMENTS/FRAGSTATS_DO CUMENTS.HTML HTTPS://WWW.R-PROJECT.ORG/OTHER-DOCS.HTML https://simlander.wordpress.com/ https://www.unocero.com/ciencia/de-fractales-y-automatas-celulares/ https://www.xatakaciencia.com/matematicas/que-son-los-fractales-y-como-se-construyen https://www.xatakaciencia.com/matematicas/que-son-los-fractales-y-como-se-construyen https://www.ibiblio.org/lifepatterns/october1970.html https://es.wikipedia.org/wiki/Juego_de_la_vida https://www.sciencedirect.com/science/article/pii/S0198971519300948?via%3Dihub https://www.sciencedirect.com/science/article/pii/S0198971519300948?via%3Dihub 46 8- ANEXOS Evolución de suelo urbano año 2000-2009 sobre la CAM. CALIBRACIONES DEL AC EN R Matrices utilizadas en vecindad w (0,0,50,0,0, 0,50,50,50,0, 50,50,500,50,50, 0,50,50,50,0, 0,0,50,0,0) w2 (0,0,5,0,0, 0,5,5,5,0, 5,5,50,5,5, 0,5,5,5,0, 0,0,5,0,0) w3(0.000,0.000,0.500,0.000,0.000, 0.000,29.496,50.000,29.496,0.000, 0.500,50.000,500.000,50.000,0.500, 0.000,29.496,50.000,29.496,0.000, 0.000,0.000,0.500,0.000,0.000), Cálculos usados en accesibilidad Ac1 (1+(distroad/1))^-1 Ac2 (1+(distroad/1000))^-1 Ac3 (1+(distroad/1000000))^-1 El orden de las calibraciones basadas en la vecindad es de : W1-W2-W3. 47 Accesibilidad 1 48 Accesibilidad 2 49 Accesibilidad 3 50 MAP COMPARISON KIT Mapa de observación de diferencias en los usos del suelo comparados entre 2000- 2009 y 2009 simulado con AC. En estos mapas se observan las diferencias del mapa simulado en comparación con los años 2000 y 2009. Accesibilidad 1 → W1 → W2 → W3 Accesibilidad 2 → W1 → W2 → W3 51 Accesibilidad 3 → W1 → W2 → W3 Accesibilidad 4 → W1 → W2 → W3 52 Script desarrollado. Ejemplo de AC con Accesibilidad 1 y vecindad 2; modelo estático. #begin script #set workspace and libraries setwd("C:/simlandertfm/SIMLANDERICHARD") require(raster) set.seed(102018) library(sf) library(igraph) #begin import data mlu00 <- raster("mlu_20001.tif") print("imported mlu00") mlu00[mlu00%in%c(1)] <- 0 mlu00[mlu00%in%c(11)] <- NA mlu00[mlu00 > 0] <- 1 lu1 <- mlu00 mlu09 <- raster("mlu_20091.tif") print("imported mlu09") mlu09[mlu09%in%c(1)] <- 0 mlu09[mlu09%in%c(11)] <- NA mlu09[mlu09 > 0] <- 1 lu2 <- mlu09 slope <- raster("recslope_cam.tif") recslope <- reclassify(slope,c(-Inf,1,1, 1,2,0.9, 2,3,0.8, 3,4,0.5, 4,Inf,0.1)) recslope <- resample(recslope, mlu00, method ="bilinear") #suitability <- (recslope) distroad <- raster("4_dist_roads100.asc") distroad <- resample(distroad, mlu00, method ="bilinear") #end import data #----------------------- #begin tidy up a bit - calculate a change map, set null values for features and zoned areas distroad_null <- mask(distroad,lu1) slope_null <- mask(recslope,lu1) #----------------------- #begin calculate accessibility. #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ ###ACCESIBILIDAD SEGUN WHITE ET AL 1997, DISTINTOS VALORES DE COEFICIENTE accurb <- (1+(distroad/1))^-1 accurb2 <- (1+(distroad/1000))^-1 accurb3 <- (1+(distroad/1000000))^-1 #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ 53 ###OTRA VARIANTE PARA CALCULAR ACCESSIBILITY SEGUN ROODPOSHTI ET AL 2020 (VER A https://www.sciencedirect.com/science/article/pii/S0198971519300948) #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ #begin calculate accessibility ################################ # !!!!!!!!!!!!!!!!!!!!!!!!!!!! # ################################ # The crisp version is not smooth enough # This sigmoid function will help to achieve a less cluttered simulation con <- function(condition, trueValue, falseValue) { return(condition * trueValue + (!condition)*falseValue) } # this is sigmoid function using con # for scoring accessibility accurb4 <- con( (distroad < 300), 1, con( (300 <= distroad) & (distroad < 1500), (0.5*(1 + cos((3.14159*(distroad-300))/1200))),0)) # end calculate accessibility #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ model_accessibility <- accurb #we have a normalised accessibility map #acc<-calc(accessibility, funacc) #end calculate accessibility #----------------------- #begin calculate suitability #slope_suit <- reclassify(slope_null,c(- Inf,0,1,0,5,0.9,5,10,0.8,10,15,0.5,15,20,0.3,20,Inf,0.1)) model_suitability <- recslope #end calculate suitability #----------------------- #set weights matrix for neighbourhood #w <- matrix(c(0,0,50,0,0, # 0,50,50,50,0, # 50,50,500,50,50, # 0,50,50,50,0, # 0,0,50,0,0), # w <- matrix (c(0.000,0.000,0.005,0.000,0.000, 0.000,0.031,0.050,0.031,0.000, 0.005,0.050,0.500,0.050,0.005, 0.000,0.031,0.050,0.031,0.000, 0.000,0.000,0.005,0.000,0.000), nr=5,nc=5) w2 <- matrix(c(0.000,0.000,0.500,0.000,0.000, 0.000,3.136,5.000,3.136,0.000, 0.500,5.000,50.000,5.000,0.500, 0.000,3.136,5.000,3.136,0.000, 0.000,0.000,0.500,0.000,0.000), nr=5,nc=5) w3 <- matrix(c(0.000,0.000,0.500,0.000,0.000, 0.000,29.496,50.000,29.496,0.000, 54 0.500,50.000,500.000,50.000,0.500, 0.000,29.496,50.000,29.496,0.000, 0.000,0.000,0.500,0.000,0.000), nr=5,nc=5) #----------------------- #preparations for simulation Tn <- lu1 #set T1 to first map lu1. Change this if we want to start from a different map. print (paste("simulation starting at ", Sys.time(), sep="")) startyear <- 2000 endyear <- 2009 total <- (endyear)-(startyear) #-------------------------------calculate demand from lu1 and lu2 (must have 2 rows and 2 columns with urban land in row 2, column 2) dflu1 <- as.data.frame(freq(lu1)) dflu2 <- as.data.frame(freq(lu2)) urbdemand <- as.numeric(dflu1[2,2]) finaldemand <- as.numeric(dflu2[2,2]) andem <- (finaldemand - urbdemand)/total andem <- round(andem,0) #-------------------------------end calculate demand from lu1 and lu2 print("set demand, starting timer..") ptm <- proc.time() #end preparations #----------------------- pb <- txtProgressBar(min = 0, max = total, style = 3) for (i in 1:total){ #begin running simulation urbdemand1 <- urbdemand + andem*i #begin neighbourhood block n <- focal(Tn, w=w2, na.rm=TRUE) #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ "na.rm=TRUE" MUY IMPORTANTE PARA NO TENER EFECTO "CUADRADO" nhood <- cover(n, Tn) model_nhood <- nhood #end neighbourhood block #----------------------- #begin random block x <-runif(ncell(Tn)) #corrected to calculate directly from number of cells in the map weibull <- 1+(-log(1-(x)))*exp(1/2) funselect <- function(x) { x[x!=0] <- NA; return(x) } #extract only vacant areas (value 0) so that existing functional land uses are not randomized. Set everything else to NA vacants <- calc(Tn, funselect) #BE CAREFUL! NEED TO TAKE VACANTS FROM T1, NOT LU56 random <- lu1 #just copying lu1 as a skeleton map values(random)<-weibull model_random <- mask(random, vacants) #generating a mask to apply NAs to the random layer. model_random <- cover(model_random,Tn) #cover to fill the NAs back in with the original values from T1. #end random block #----------------------- #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ #begin dynamic accessibility block 55 #print("calculating dynamic accessibility.......") #model_accessibility <- accurb/i ##AQUI COGEMOS EL MAPA DE ACCESSIBILIDAD ORIGINAL Y LO DEBILITAMOS SEGUN AVANZA LA SIMULACIÓN #d <- gridDistance(Tn, origin = 1) ###TIENES QUE INSTALAR igraph - install.packages("igraph") #d2 <- abs(d-(cellStats(d,max))) ##INVERTIR efecto de distancia para dar valores mas altos a zonas mas cercanas al uso urbano #d3 <- d2/(cellStats(d2,max)) #normalizar #d4 <- model_accessibility + (d3/1000) #RECALCULAR ACESSIBILIDAD #d5 <- d4/(cellStats(d4,max)) ##NORMALIZAR #model_accessibility <- d5 #@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ #end dynamic accessibility block #begin transition potential calculation model_TTP <- (model_accessibility+1)*(model_suitability+1)*(model_nhood+1)*(model_r andom) fun <- function(x) { x[is.na(x)] <- 0; return(x)} TTP <- calc(model_TTP, fun) #Select the highest 164 values from the TTP map #r1demand <- 164 #this is the number of cells we want to allocate from r1 x <- as.matrix(TTP) n <- urbdemand1 #the demand calculated as above x2 <- sort(-x, partial = n) x2h <- -sort(x2[1:n]) ix <- which(x %in% x2h) rowsT1 <- nrow(Tn) colsT1 <- ncol(Tn) ro <- ix %% rowsT1 ro <- ifelse(ro == 0L, rowsT1, ro) co <- 1 + ix %/% rowsT1 x3 <- x[ix] d <- data.frame(row = ro, col = co, x = x3) result <- d[rev(order(d$x)), ] #----------------------- #test for duplicates that inflate the number of cells allocated difftrans <- (length(result$x)-n) if (difftrans > 0) { result2 <- head(result,-difftrans) #remove the duplicates from the end of the file (the weakest candidate cells) result <- result2 } #turn selected n values in the dataframe into a matrix and then to a raster x.mat <- matrix(0, rowsT1, colsT1) #create a matrix with the right number of rows and columns and fill with 0 values #x.mat[cbind(result$row,result$col)] <-result$x #works just fine. x.mat[cbind(result$row,result$col)] <-1 #but actually we want values to be 1 (urban land), not the TP value. #which(!is.na(x.mat), arr.ind =TRUE) #pick out the non-NA values. Useful for testing that this has actually worked. r <- raster(x.mat) extent(r) <- extent(Tn) 56 newdata <- r newdata <- mask(newdata, lu1) #generating a mask to apply NAs to the final map layer. Tn <- newdata #directly allocated all the cells at their most favourable locations according to TTP filen <- paste("lu", (startyear + i), ".png", sep="") #filenasc <- paste("lu", (startyear + i), ".asc", sep="") #LUout <- writeRaster(T1, filename=(filenasc), format="ascii", overwrite=TRUE) #filen <- paste(filen,".png", sep="") png(filename = paste(filen),width = 1200, height = 1200, bg="white") #Plot each land use map to be able to make an animation. #plot(T1) #plot(T1,breaks=breakpoints,col=colors,main=(startyear+ii)) plot(Tn,main=(startyear+i)) dev.off() removeTmpFiles(h=5) Sys.sleep(0.1) setTxtProgressBar(pb, i) print (paste("FINISHED: landuse simulation for ", (startyear + i), sep="")) print(proc.time() - ptm) } close(pb) png(filename = paste("lu",startyear,".png",sep=""),width = 1200, height = 1200, bg="white") #Plot initial map also to animations file plot(lu1, main = paste("initial land use map - ",startyear,sep="")) #lu1 for calibration, lu2 for simulation phase dev.off() #eval_anim <- paste("system ('convert -delay 200 -loop 0 *.png lu",startyear,"_",endyear,".gif')",sep="") #creates an animated gif using imagemagick, comment back in if you are on a linux system and have imagemagick installed #eval(parse(text=eval_anim)) #comment back in if you are on a linux system and have imagemagick installed #simlu2011 <- T1 sim09 <- Tn writeRaster(sim09,"sim09_ac1_w2_smooth.asc",format="ascii",overwrite=T ) sss <- stack(lu2,sim09) #name1 <- paste("initial land use ",startyear,sep="") name1 <- paste("real land use ",endyear,sep="") name2 <- paste("simulated land use ",endyear,sep="") names(sss) <- c(name1,name2) print(plot(sss)) #without the print command, the output will not appear on the screen print (paste("simulation complete at ", Sys.time(), sep="")) #print(plot(sss,breaks=breakpoints,col=colors)) #same as above, but just in case you have colors and breakpoints defined