Introducción

La alteración de la superficie terrestre, a causa de las acciones humanas, ha generado que cambien los usos del suelo en la Tierra. Por esta razón, la caracterización de las coberturas y usos del suelo es clave para comprender las interacciones antrópico-ambientales a diferentes escalas espacio-temporales (^{Müller & Munroe, 2014}). Los patrones espaciales y sus cambios tienen implicaciones significativas cuando se analizan, por ejemplo, el balance hídrico, el cambio climático, la deforestación, la sedimentación, la pérdida de biodiversidad, y los procesos biogeoquímicos de la Tierra (^{DeFries & Eshleman, 2004}).

En las últimas décadas se ha observado el aumento acelerado de la población mundial, acompañado de la intensificación de actividades agrícolas y de expansión de la frontera agropecuaria (^{Eitzinger, Orlandini, Stefanski, & Naylor, 2010}). Estas transformaciones estructurales han provocado la alteración de procesos hidrológicos (^{Breuer et al., 2009}; ^{Fohrer, Haverkamp, Eckhardt, & Frede, 2001}; ^{Li, Liu, Zhang, & Zheng, 2009}), desbalances en los ciclos del carbono y nitrógeno (^{Murty, Kirschbaum, Mcmurtrie, & Mcgilvray, 2002}; ^{Post & Kwon, 2000}) y pérdidas de especies y ecosistemas a tasas altas, a pesar de los esfuerzos por conservar la diversidad biológica (^{Pettorelli et al., 2016}).

En Argentina, particularmente en la región pampeana, se ha experimentado un proceso de agriculturización caracterizado por el incremento de las áreas agrícolas en zonas de uso tradicionalmente ganadero. Dicho proceso comenzó en la década de 1960 y se ha acelerado de manera dramática durante los últimos 20 a 25 años, sostenido principalmente por cambios económicos estructurales y la incorporación del paquete tecnológico asociado al cultivo de soya (soja en Argentina) genéticamente modificada (^{Manuel-Navarrete et al., 2009}). Estos cambios impactan fuertemente la estructura y funcionamiento de los ecosistemas, se perciben a diferentes escalas (^{Baldi, Guerschman, & Paruelo, 2006}; ^{Guerschman, Paruelo, & Burke, 2003}; ^{Lara & Gandini, 2014}) y generan alteraciones en los patrones de productividad, estabilidad y sostenibilidad (^{Viglizzo, Roberto, Lértora, Gay, & Bernardos, 1997}).

La disponibilidad y uso de imágenes satelitales han permitido el análisis de grandes extensiones de superficie, adquiriendo información a escalas temporal y espacial muy diferentes (^{Horning, Robinson, Sterling, Turner, & Spector, 2010}). El uso de dichas imágenes está asociado a la descripción de tipos de coberturas/usos del suelo, mediante la caracterización de la heterogeneidad del paisaje y su patrón espacial (^{Cabello & Paruelo, 2008}; ^{Paruelo, 2008}). Tal aplicación está sustentada en las formas de respuesta de los elementos de la superficie terrestre a distintas longitudes de onda del espectro electromagnético (^{Chuvieco, 1990}), aunque a menudo una de las limitantes principales es la escasa información de campo o su captura inadecuada. En este contexto, varios índices espectrales derivados de sensores remotos han sido propuestos como estimadores de la presencia y condición de la vegetación, basados en la observación de las firmas espectrales de reflectancia típica de las hojas. En estas, la energía reflejada en el rango visible del espectro es muy baja, debido a la alta absorción de los pigmentos fotosintéticamente activos, con absorciones máximas en las longitudes de onda del azul (470 nm) y del rojo (670 nm); casi la totalidad de la radiación del infrarrojo cercano (IRC) es dispersada (reflejada y transmitida) con muy poca absorción, debido a la estructura del mesófilo de las hojas (^{Solano, Didan, Jacobson, & Huete, 2010}). Como resultado, el contraste entre las respuestas del rojo y del IRC es una medida sensible de la cantidad de vegetación, con máximas diferencias sobre la vegetación densa y mínimas donde existe poca o nula cobertura (Solano et al., 2010). Uno de los índices espectrales más utilizados es el índice de vegetación de diferencia normalizada (NDVI, por sus siglas en inglés), el cual es un estimador lineal de la fracción de radiación fotosintéticamente activa interceptada por la vegetación (fAPAR) (^{Wang et al., 2004}). A partir de los datos brindados por los sensores satelitales de alta resolución temporal (tiempo corto de revisita, i. e. MODIS) pueden construirse series hipertemporales de NDVI, de las cuales es posible extraer información valiosa acerca de la fenología de la vegetación (^{Reed, Schwartz, & Xiangming, 2009}). De esta manera, es posible utilizar las diferencias en los patrones estacionales para caracterizar diversos tipos de vegetación (cultivos de verano, invierno, vegetación perenne, etc.). Por otro lado, los mapeos disponibles de coberturas y usos del suelo a gran escala (^{Blanco et al., 2013}; ^{Eva et al., 2004}) presentan inconsistencias debido a su aproximación continental y resultan insuficientes para poder ser aplicados a una escala más detallada. En este sentido, el objetivo del trabajo fue elaborar una metodología para la caracterización del uso del suelo en la cuenca superior del arroyo del Azul en Buenos Aires, Argentina, mediante el uso de imágenes de media resolución espacial de los satélites Landsat y MODIS.

Materiales y métodos

Área de estudio

La cuenca superior del arroyo del Azul (Figura 1) se ubica en el centro de la provincia de Buenos Aires, Argentina, entre 59° 48’ - 60° 0’ LO y 37° 01’ 48’’ - 37° 09’ 36’’ LS, drenando una superficie de 1 024 km². Dicha cuenca nace hacia el Sur en la localidad de Chillar (60 km al sur de la ciudad de Azul); sus principales afluentes son el arroyo Videla con una superficie de 120 km² y el de Santa Catalina con 140 km².

Figura 1 Cuenca superior del arroyo del Azul ubicada en el centro de la provincia de Buenos Aires, Argentina. 

La altitud en la cuenca varía de 142 a 367 m. La parte más alta, ubicada hacia el sureste, se encuentra sobre el Sistema de Tandilia, donde las pendientes son mayores de 6 %. Hacia el noroeste, la cuenca se transforma en una llanura de transición, en el sector medio-bajo, con pendientes menores de 1 %. La cuenca limita al Norte con la ciudad de Azul, al Oeste con la ciudad de Olavarría, al Este con la ciudad de Tandil y al Sur con la ciudad de Chillar. La precipitación media anual es de 913 mm. Las mayores precipitaciones se presentan en el periodo primavera-verano con un promedio mensual alrededor de 100 mm, y las más bajas ocurren en otoño-invierno con una media mensual de 50 mm. Los meses de diciembre a marzo son los de mayor temperatura con un promedio de 20 °C, mientras que los meses de junio a agosto son los más fríos con un promedio de 8 °C.

Clasificación y obtención de coberturas del suelo

A partir de la serie temporal de imágenes NDVI-MODIS se hizo una clasificación no supervisada con el algoritmo ISODATA (Iterative Self-Organizing Data Analysis), el cual hace la partición de una imagen de n-dimensiones en un número de conglomerados de acuerdo con los parámetros establecidos. Por lo general, el proceso comienza dando un número de centroides aleatorios y, posteriormente, cada uno de los pixeles es asignado a aquel conglomerado cuyo centroide se encuentra más cercano, en términos de las distancias euclidianas. Este es un proceso iterativo en el cual las medias (centroides) se recalculan y los pixeles se reclasifican con respecto a las nuevas medias. El proceso continúa hasta que el número de pixeles en cada clase cambia menos que el umbral de cambio propuesto o hasta que se alcanza el número máximo de iteraciones configurado (^{Tou & Gonzalez, 1974}).

En este procedimiento se consideraron 50 iteraciones, un umbral de tolerancia del 5 % y desviación estándar máxima de 1. De esta forma, se generaron 15 clases de las cuales se obtuvo el valor medio de NDVI para cada una de las 23 fechas contempladas en la serie temporal. Los valores medios corresponden a las firmas fenológicas, las cuales fueron asociadas, con base en el comportamiento anual, a cada una de las coberturas y usos del suelo definidas: soya, maíz, sistemas de doble cultivo trigo-soya, pasturas (naturales y cultivadas) y áreas mezcladas (Figura 2). Por ejemplo, los cultivos se asociaron a las clases que presentaron uno o dos períodos de marcado crecimiento y un período con bajos valores de NDVI. Para discriminar entre cultivos, las curvas unimodales con valores máximos de NDVI en febrero se clasificaron como soya, y las que mostraban un patrón similar, pero con los valores máximos de NDVI en enero y febrero se clasificaron como maíz. Aquellas clases que presentaron un máximo durante septiembre-octubre (característico de cultivos invernales) y, además, un máximo en febrero, se categorizaron como trigo-soya.

Figura 2 Curvas estacionales del índice de vegetación de diferencia normalizada (NDVI) en el período mayo del 2015 a mayo del 2016 para cada una de las 15 clases resultantes de la clasificación no supervisada, escena h13v12 del satélite MODIS. a) área mezclada, b) soya, c) maíz, d) trigo-soya y e) pastura.  

Por otro lado, las clases que no presentaron un patrón promedio unimodal o bimodal se clasificaron como pasturas y áreas mezcladas. Las pasturas, que incluyen pastizales perennes cultivados y pastizales naturales, se caracterizaron por mantener un valor promedio de NDVI mayor de 0.4 a lo largo de la estación de crecimiento; mientras que las áreas mezcladas corresponden a un valor estacional promedio menor de 0.4. Estas últimas corresponden a zonas urbanas, cuerpos de agua y áreas de afloramiento rocoso con vegetación activa.

Es probable que las áreas mezcladas que no lograron diferenciarse a partir de las imágenes MODIS estén asociadas a pixeles mixtos a esa escala, lo cual dificultaría hallar pixeles puros para diferenciar entre coberturas (^{Lara & Gandini, 2016}). Por esta razón, para poder discriminar las áreas, se utilizaron dos imágenes Landsat 8 (época invernal y estival) en una combinación 7, 6, 4, que permite penetrar partículas atmosféricas, humo y neblina (Figura 3). Posteriormente, las dos fuentes de datos se fusionaron remuestreando la clasificación realizada con imágenes MODIS previamente, a una resolución espacial de 30 m, para obtener el mapa final de coberturas y usos del suelo de la cuenca superior del arroyo del Azul.

Figura 3 Imágenes Landsat 8 utilizadas. (a) Época invernal, 15 de agosto de 2015. (b) Época estival, 22 de enero de 2016. Path/Row 225/086. Composición de bandas RGB: 7, 6, 4. La vegetación aparece en tonos de verde oscuro; las áreas urbanas son de color gris, azul verdoso o púrpura; el agua es de color negro o azul oscuro; y los afloramientos rocosos aparecen en tonos de rojo o amarillo. 

Resultados y discusión

La Figura 4 muestra el mapa final de coberturas y usos del suelo de la cuenca superior del arroyo del Azul, Argentina. El mapa se obtuvo a partir de la combinación de imágenes de media resolución espacial, para discriminar tipos de cultivos y áreas de pastizales con base en sus patrones fenológicos característicos. El Cuadro 1 indica la superficie cubierta por cada uno de los usos del suelo; el sistema de doble cultivo trigo-soya es el más extendido representando 39.4 % (40 470.4 ha) del área de estudio. La suma de las superficies cubiertas por todos los usos agrícolas alcanza 81.5 %, lo cual refleja la magnitud del proceso de agriculturización en la región, teniendo en cuenta que la vegetación natural era el pastizal natural (^{Matteucci, 2012}). A pesar de que no existen datos oficiales, las estadísticas agropecuarias del Ministerio de Agroindustria de Argentina (^{Sistema Integrado de Información Agropecuaria [SIIA], 2016}) reflejan en cierta medida, el predominio de los mismos cultivos en la campaña 2015-2016 para el departamento de Azul, tales como soya (124 700 ha), trigo (132 110 ha) y maíz (43 300 ha); se debe tener en cuenta que en estos datos existen áreas de uso agrícola que son externas al área de estudio.

Figura 4 Mapa final de usos del suelo en la cuenca superior del arroyo del Azul (Buenos Aires, Argentina) para el periodo 2015-2016, mediante el uso de imágenes MODIS y Landsat 8.  

En general, las áreas cultivadas se encuentran en la porción más alta de la cuenca, ya que ahí predominan los suelos de mayor aptitud y con pocas restricciones para la actividad agrícola, aunque en los últimos años, el avance de los cultivos también se ha dado sobre áreas de uso tradicionalmente ganadero (^{Lara & Gandini, 2014}). Las áreas de pasturas (que incluyen a los pastizales naturales y las pasturas cultivadas) se encuentran asociadas, principalmente, a los afloramientos rocosos del Sistema de Tandilia y en zonas bajas de la cuenca. Probablemente, las zonas cercanas a los afloramientos rocosos representan zonas de pastizal natural donde no es posible labrar el suelo, mientras que las ubicadas en zonas bajas de la cuenca representan, en su mayoría, pasturas cultivadas que suelen utilizarse en zonas con limitaciones de drenaje y salinización. La metodología actual no permitió obtener una separación clara en las firmas fenológicas de estos pastizales, debido a que el solapamiento de las fechas de inicio y final de las estaciones de crecimiento fue frecuente.

El mapa de usos del suelo mostró buenos niveles de aciertos, con una precisión total de 88.9 % y un coeficiente Kappa de 0.86. De acuerdo con el Cuadro 2, las precisiones más bajas del productor correspondieron a los usos maíz y soya con valores de 69.9 y 76.0 %, respectivamente. Justamente, la mayor confusión en la clasificación se da entre estos dos usos, lo cual es esperable ya que sus firmas fenológicas son similares (Figura 2). A pesar de esto, la precisión del mapa final, obtenido a partir de una metodología basada en firmas fenológicas derivadas de NDVI y sin información de campo, representó de manera satisfactoria los patrones espacio-temporales de los principales usos de la cuenca en un período determinado. Sin duda, la combinación de datos obtenidos a una mayor resolución espacial contribuyó a lograr una precisión mayor.

La metodología desarrollada posee la ventaja de ser relativamente rápida y replicable con pocos insumos, para ser utilizada en diferentes regiones con el fin de mapear los usos de mayor representatividad; sin embargo, tiene la desventaja de que el comportamiento anual de los usos a discriminar se debe conocer a priori, lo que podría evitarse tomando datos de campo en pixeles puros para obtener la firma fenológica.

La caracterización de las coberturas y usos del suelo representa uno de los insumos imprescindibles para el monitoreo y manejo de los recursos naturales a diferentes escalas, como por ejemplo en la modelación hidrológica semidistribuida (^{Arnold et al., 2012}); sin embargo, en Argentina suele ser una de las informaciones ausentes en los organismos de gestión departamental o regional. Por lo tanto, cuando se requieren estrategias de manejo sustentable de los recursos naturales a escala local, los procedimientos existentes a gran escala (^{Blanco et al., 2013}; ^{Eva et al., 2004}) resultan deficientes, y es entonces cuando la metodología desarrollada en este trabajo puede aportar herramientas útiles. Por otro lado, en regiones de gran dinamismo de usos del suelo, como la aquí analizada, la metodología permitiría evaluar y monitorear los cambios de coberturas y, en particular, de la frontera agropecuaria, prácticamente en tiempo real.

Conclusiones

La metodología desarrollada en el trabajo permitió la caracterización de las coberturas y usos del suelo de la cuenca superior del arroyo del Azul, Argentina, mediante el uso de firmas fenológicas derivadas de series temporales de NDVI-MODIS. Las áreas mezcladas que no pudieron diferenciarse fueron categorizadas a partir de imágenes Landsat 8, para mejorar la precisión del mapa final. El método de clasificación empleado es una alternativa válida cuando no se dispone de datos precisos de campo, pero requiere conocimiento acerca del comportamiento estacional de las coberturas a ser mapeadas. De esta manera, la propuesta presentada puede ser aplicable y reproducible en diferentes regiones y escalas, convirtiéndose en una herramienta valiosa para el monitoreo de los cambios de usos del suelo a escala local y regional.