AUTOMATED MAPPING OF PERMANENT PRESERVATION AREAS ON HILLTOPS

Permanent Preservation Areas (PPAs) on hilltops are among the many areas protected by the New Forest Code in Brazil. Mapping of these involves difficult interpretation and application of the Law, as well a complex task of translating it in map algebra. This paper aims to present, in detail, a methodological model for delimitation of PPAs on hilltops, according to the Brazilian new Forest Code (NFC, Law 12,651/2012). The model was developed in Model Builder for ArcGIS 10.2, and is able to map the PPAs in any digital terrain model. However, field validations are required to verify its efficiency. There is need for legal standardization of criteria that may cause subjectivity in delimitation. The organization of these data on a large scale is very important, as example, to the Rural Environmental Registry, which provides georeferencing of all rural properties and its protected areas in Brazil. MAPEAMENTO AUTOMATIZADO DE ÁREAS DE PRESERVAÇÃO PERMANENTE EM TOPO DE MORROS RESUMO: As Áreas de Preservação Permanente em topos de morros estão entre as diversas áeras protegidas pelo novo Código Florestal (NCF) brasileiro. O mapeamento das mesmas envolve dificuldades técnicas de interpretação da Lei, bem como a complexa tarefa de traduzi-la em álgebra de mapas. O presente trabalho apresenta, em detalhe, um modelo metodológico para delimitação de APPs em topos de morros, de acordo com o NCF (Lei 12.651/2012). O modelo foi construído em ambiente Model Builder para ArcGIS 10.2 e é capaz de mapear as APPs em qualquer modelo digital de terreno. No entanto, validações em campo são necessárias para verificar a eficiência do modelo. Também, é necessária a regulamentação de critérios que causam subjetividade na delimitação. A organização desses dados em larga escala é importante, por exemplo, para o Cadastro Ambiental Rural (CAR), que prevê o georreferenciamento de todas as propriedades rurais e suas áreas protegidas, no Brasil.


INTRODUCTION
Remote sensing and geographic information systems are essential tools for environmental management and monitoring natural resources.As advantages, it brings greater speed to produce georeferenced data and maps.The use of geoprocessing techniques for delimiting protected areas contributes to generation of geo-referenced databases, which are important information providers to the state research institutions and the private sector.
In Brazil, the management of natural vegetation in private properties is guided by a set of rules known as the Forest Code.The first version of the code is from 1965 (Law 4,771), which was largely amended over the subsequent decades via decrees and regulations.The New Forest Code (NFC), Law 12,651/2012(BRASIL, 2012), is an updated and compiled version of the previous legislation (SOARES-FILHO et al., 2014).
The NFC provides the status of Permanent Preservation Area (PPA) to certain areas, due to their importance in protecting soil and water, as well as the most unstable places of diverse landscapes of the country.In these places, the natural vegetation should be kept to mitigate soil erosion, provide greater water production and maintenance of its quality, among other environmental benefits like protecting biodiversity and sequestrating carbon dioxide.Riparian zones, water springs, steep slopes and mangroves are examples of PPAs.
The mapping of PPAs integrates a set of core analysis for planning and management of natural resources and territory (PEDRON et al., 2006;SOARES et al., 2011).The Rural Environmental Registry (CAR, acronym in Portuguese) established by the NFC, aims to create a database of land use and natural sites of all farms in the country.This remarkable challenge will be achieved, largely, thanks to geographic information systems.
Unlike other PPAs, the delimitation of the hilltops is quite complex (RIBEIRO et al., 2005), both in the field and in the laboratory.Consequently, PPAs on hilltops are often set aside in environmental spatial analysis.Even in the CAR, it is expected an underestimation of protected areas on hilltops, since many will not be declared by the users because of the difficulty of application of the rule.In addition, the Government still does not have a tool for verifying undeclared areas.
The high financial cost and the time demand for the on-site mapping justify the development of automated digital mapping techniques.The automated process eliminates subjectivity, while suppressing the interference of the analyst in situations that could cause dubiety, since these decisions will be made by the algorithm.Thus, mapping several times in the same place, using the same basic topographic data, will always return the same output.
The mapping of PPAs on hilltops is based on digital terrain models and geoprocessing tools (RIBEIRO et al., 2005;SOARES et al., 2011;VICTORIA et al., 2008).Although digital representations of relief are approximations of reality, and are subject to errors and uncertainties (FISHER and TATE, 2006), the development of generic models of automated mapping represent major advance on the issue.
Based on the above, the objective of this study was to develop a new methodology for automated mapping of PPAs on hilltops in wavy relief, according to the new Brazilian Forest Code.The methodological model was built from the legal requirements translated in a series of procedures that uses digital terrain models as input data.

METHODOLOGY
The model was developed in Model Builder programming environment, ArcGIS 10.2 system, which operates in Python language.To serve the growing demand of free software use, operations were explained in order to enable its replication in other geographic information systems, using similar functions.
All functions were performed following the provisions of Article 4, section IX of Law 12,651 / 2012(BRASIL, 2012), for hilltops in wavy relief.Mapping of hilltops in plan relief is not addressed in this work.The Article 4 says (in translation): "IX -on top of hills and mountains, with a minimum height of one hundred (100) meters and an average slope greater than 25°, the areas bounded from the corresponding contour of 2/3 (two thirds) of minimum elevation always in relation to the base, which is defined by the horizontal plane determined by adjacent plain or water surface or, in wavy relief, the height of the nearest saddle point" (BRASIL, 2012).
The model was built in a modular structure, dividing the geoprocessing operations in groups with independent functions, but with a complete sense within each one (e.g., a part of the model defines the bases of the hills, the other identifies the tops, and so on).The processing was divided into nine different steps, where the products generated in each step are used as input data in the subsequent steps.The operations performed in GIS were organized in tables with the purpose of systematizing and ease the replication of the methodology.
As input data, any Digital Elevation Model (DEM) grid in UTM projection can be used with this methodology.The Law 12,651 / 2012 does not define, so far, a scale for mapping PPAs on hilltops, nor which database should be officially adopted.Althought theoretically any DEM could be used, Victoria (2010) showed that the scale of the cartographic data affected the result of the mapping, considering the earlier legislation (Law 4,771) to map the PPAs.
To expose the results of the process, an ASTER DEM with 30-meter spatial resolution was used and a study area with strong relief was selected.This area is a portion of the Serra do Gandarela, located in Minas Gerais State, Brazil, with an extention of 40 x 40 km.The final output of the process is a vector file (shapefile), in which each feature is a PPA in hilltop.

RESULTS AND DISCUSSION
A model has been successfully built in Model Builder environment.The outputs of the processes are consistent with the rules defined by the law and the processing time is feasible.In all operations where the setting were applicable, the Environment Options were enabled to limit the extent of processing to the DEM, as well as adherence (snapping) to the cells of this DEM.Thus, all files of the matrix type (raster) that have been created, had their cells aligned with the original DEM.

Separating the hills by a hydrological concept
The first step was to map the areas that define the extent of each hill (Table 1).Thus, the hills are divided, considering that its limits are defined by the saddle points, and also matched with the drainage network.This area is based on the contribution area relating to inverted relief.This limit is not confused with the one defined by law as hill (whose base is defined by the altitude of the closest saddle point of a hill), nor with the Permanent Preservation Area.To get to the delimitation of the area covered by the law as hill, various functions were carried out, based on the delimitation of what might be called a simplistic hydrological base of the hills, which are coincident with hydrography and saddle points.For a better understanding, the analyst can overlay the layer produced at the end of this section with contour lines, where one will notice the coincidence of the lines that divide the hills with saddle points.
In some situations, there is a need of a brief data pre-processing before the execution of the functions, particularly in DEMs presenting a strong microrelief expression, such as the high resolution ones (CHAPLOT et al., 2006).In these instances, the need for pretreatment comes from the fact that tiny relief imperfections could be considered as saddle points, which has no practical sense.To work around this problem, an averaging filter in rectangular moving window of 3x3 size was used (FOCAL STATISTICS [1.1]).Then spurious depressions (sinks) were removed with the FILL command [1.2].The use of filtering techniques of the terrain artifacts is common (ALBANI and KLINKENBERG, 2003) and provides creation of more consistent products by eliminating systematic errors that lead to abrupt variations in the data locally.Producing matrices with the average values of a cell and its neighborhood is an underlying trend of the terrain parameters (CARLISLE, 2005).Although the pre-treatment is necessary, the DEM produced by it will be used only in the first step of the methodology, where a previous separation of the hills is made.All altitude values used in the model come from the original DEM.
Using the pre-processed DEM, the next step aims to numerically invert the relief, with a subtraction operation applied by MINUS [1.3] command.In the generated raster, the value of each cell will be equal to "8,000" minus the value indicated in the previous DEM for that cell.Subsequently, the flow direction of this inverted DEM was generated (FLOW DIRECTION [1.4]), which was further used in BASIN command [1.5] to delimit the drainage basins, which are delineated by the ridge line therebetween.These lines, which pass through the inverted DEM top points, when overlapping the normal relief, pass exactly over the saddle points.The matrix file produced by BASIN command was then converted to type integer (INT [1.6]) and finally converted into a polygon (RASTER TO POLYGON [1.7]).The final aspect of this section's product is illustrated in Figure 1.

Finding the tops of the hills
Once the hills were individualized, the next step was to identify the top point in each of them.Thus, for each area delimited in the previous section, only a single top point defined by the maximum altitude present in each of these areas was identified.Besides the original DEM, the output file generated in the function RASTER TO POLYGON [1.7], which contains the boundary of each of the hills, was used.
The first function performed in this step (ZONAL STATISTICS [2.1]), was to calculate, based on the edge of the hills, the maximum altitude value (represented in DEM) within these areas, where the maximum value was assigned to all cells in the zone, in a new raster file ("basin_max", as seen in Table 2).Then, the raster file was equaled to DEM by RASTER CALCULATOR [2.2], so that the cells that assume true values for expression were marked with the value 1, and the others with value 0. Within an area that defines the hill, there will be a single point where the maximum altitude value in that area is equal to the altitude of the MDE, which are precisely the top point.This matrix was then reclassified (Reclassify [2.3]), converting the value "0" to "NoData" and keeping the value "1" unchanged.Later, this file has been converted to the vector points format (Figure 2) with the RASTER TO POINT function [2.4] and, finally, the value of altitude was added to these points in the attribute table (ADD SURFACE INFORMATION [2.5]).

Identifying the saddle points
The saddle points were defined based on the area boundary of each hill (Figure 1), since the lines which define these boundaries pass over the saddle points (Table 3).The location of the saddle points coincide precisely with the point of maximum altitude of these lines.Once the several points of saddle, surrounding an elevation, were located, the closest one to each elevation had to be determined.
Initially, the polygons which define the limits of the hills ("basin_pol" generated in step 1.The last two steps of this section consisted of adding the respective altitude value of each saddle point to the attribute table (ADD SURFACE INFORMATION [3.7]) and then moving these saddle points to match with their respective original line, through an adherence command (SNAP [3.8]).This was done because, being generated from a raster file, these points were located in the center of cells and not above those lines.The distance used in the SNAP shall be equal to the cellsize of DEM, diagonally.Thus, points will be positioned on the lines in the shortest distance possible.Numerous saddle points were generated (Figure 3), but only those closest to each hilltop were considered.

Identification of the saddle points closest to the elevations
Once the top points of the hills and saddle points were located, we proceeded to a series of functions (Table 4), which aimed to index each top point, data from his nearest saddle point, present on attribute table 3. The closest saddle point of each elevation, will later be used to define the base of the hill, according to the New Forest Code.
First the nearest saddle point among those previously mapped was located for each elevation.Through GENERATE NEAR TABLE [4.1] function, the coordinates of the nearest saddle point of each hilltop were tabulated in a spreadsheet.These coordinates were then linked to the top points shapefile attribute table using the tool JOIN FIELD [4.2].The same tool was used then to recover the altitude value of saddle points for the table of top points.These values were obtained directly from the saddle points (JOIN FIELD [4.3]).
In short, information from saddle points was added to attribute table of top points shapefile, namely: the identifier code from the nearest saddle point, its distance, latitude, longitude and altitude.With this information gathered, it was possible to calculate the difference in altitude between the top of the elevation and its closest saddle point.
A new field in the attribute table of the top points was created (ADD FIELD [4.4]), which was used to calculate the difference in altitude between the two points (CALCULATE FIELD [4.5]).Then the information in shapefile with the top points was transferred to the polygons containing the limits of the hills (SPATIAL JOIN [4.5]) and a new geoinformation was created.
At the end of this step, three spatial information were produced: the top point, saddle points and a limit of hills containing previous information (Figure 4).

Hill basis delimitation by the altitude of the saddle point
So far, what has been called the hill limit was only a first separation of them.In this pre-defined area, all saddle points that surround an elevation or hilltop were included.However, for the purposes of the Law, the base of the hill is defined by the altitude of his nearest saddle point, which location was explained earlier.In the next step we considered, within the boundary area of each hill, only the cells where the altitude was equal or superior to that of the saddle point.
From the boundaries of the hills ("basin_join") containing the information of the top point and saddle points, a temporary file was created (MAKE FEATURE LAYER [5.1]), which allowed its conversion to the raster type (POLYGON TO RASTER [5.2]).The information used to mark the values in the new file was the altitude of the saddle point (table 5).A raster file was created, where each region, defined by the limits of each hill, had the value of the saddle point altitude in their cells.Then, the difference between the values of DEM cells and this raster was calculated (RASTER CALCULATOR FIGURE 5 Fifth step: hill base defined by the altitude of saddle point.Field Type: "Double" Field Name: "dif_H"; Expression: Join Operation: "One to one" Match Option: "Intersect" Output "saddle_prox" "point_top"* "point_top" "point_top"* "point_top"* "basin_join" [5.3]), returning value 1 only to the cells where the result was equal to or greater than zero, and zero for the others.This file, in turn, was reclassified considering only the cells with the value 1 (RECLASSIFY [5.4]) and converted to vector format with the function RASTER TO POLYGON [5.5].
The hills that share common saddle points were united as a single area.For this reason, it was necessary to intersect these areas using the limits of the hills to individualize them (INTERSECT [5.6]).Thus, a new file was created that represents the base of the hills (Figure 5), pursuant to Article 4 of Law 12,651 / 2012, set by the elevation of the nearest saddle point.

Calculating the average slope of the hills
The New Forest Code does not define the criteria to be used for calculating the average slope, which must be greater than 25 degrees for the upper third of a hill to be considered PPA.Whereas there are different methods to measure the slope, as well as different numbers of measurements may be used to obtain an average, the law opens a large space for subjectivity.
In this work, we proposed to create numerous imaginary lines connecting the top of the hill to its baseline and calculate averages of the slopes obtained from them (table 6).Thus, it ensures that the steepest slopes of a hill will always be considered in the calculation.Another way to measure the average slope would be to get the average of the pixels of the slope gradient model, present in the area understood as hill base.The following describes the processes performed to calculate the average slope of the hills and selection of those with this value greater than 25 degrees.
First, the vertices of the polygons of the hill bases were converted into points, with FEATURE VERTICES TO POINTS [6.1] tool.To these points, the XY coordinates were added in the attribute table (ADD XY COORDINATES [6.2]).In this table, a new field has been created (ADD FIELD [6.3]), where the values of the altitudes of the tops have been copied (CALCULATE FIELD [6.4]) with a new field name.This was necessary because the altitudes that were in the table have been tabulated in a field automatically created by the software.For subsequent calculations it was necessary to establish a new field and erase the previous one (DELETE FIELD [6.5]), so that subsequent calculations do not overlap the values previously calculated.
Considering that points on the hill baseline have the same altitude of the saddle point, we took the values of the raster file, which contained the altitude of saddle points and added them to the points as their value altitude (ADD SURFACE INFORMATION [6.6]).The XY coordinates of the tops were added to their attribute table (ADD XY COORDINATES [6.7]) and subsequently attached to the table of points relative to their respective baseline (JOIN FIELD [6.8]).
For each point located on the base of the hills, altitude and location information was added, both for these points as well as the respective top points.These data allowed the calculation of the tangent of the imaginary line joining the tops of the hills to the bases by applying the Pythagorean Theorem.A new field created (ADD FIELD [6.9]) was used to calculate the difference between the altitude of the top points and saddle points, and divide that value by the square root of the sum of the deviations of the distances squared (CALCULATE FIELD [6:10]).As a result, we obtained, for each point, the value given by the tangent angle defining the slope between the top and this same point, located on the hill base.Of the many points created on the baseline of the hill, only some will be sampled, as this procedure has high computational cost.
To ensure that points regarding the limits of neighboring hills were not sampled, a contracture was made in polygons that define the bases of the hills, with a distance equivalent to the size of DEM cell, using a negative distance in BUFFER [6:11] command.Then, the points located on the bases of the hills were moved to coincide with this new set limit (SNAP [6:12]), and finally, only the points that were within the polygon that define the hill base were selected, with CLIP [6:13] using a spatial tolerance of 10% of the size of DEM cell.Some points located near the drainage network, showed negative values for the tangent (which was expected because drainage has a lower height) and for this reason only those with values equal or greater than zero were selected for the calculation of slope (SELECT [6:14]).
Selected points to calculate the slope were then converted to raster format, using the values of the previously calculated tangents (POINT TO RASTER [6:15]).By the the tangent value, the angle was calculated in radians, by the function of the arc tangent.This value is in radians was multiplied by 180 and divided by π (RASTER CALCULATOR [6:16]), resulting in the slope (degrees) between each point located on the baseline and the top of the respective hill.
Then, we calculated the average of such slopes on each hill, setting the base of the hills as zone (ZONAL STATISTICS [6:17]).Finally, the raster was reclassified, exporting only the unit values (RECLASSIFY [6:18]), representing the hills which average slope is equal or greater than 25 degrees (Figure 6).deleted (RASTER CALCULATOR [7.2]) and, finally, the raster produced was reclassified (RECLASSIFY [7.3]), considering only the value "1", which represents the hills with height equal or greater than 100 meters (Figure 7).

Measuring the height of the hills
Another criteria defined by law, for a hill to have Permanent Preservation Area, is its height.Based on the altitude of the saddle point, the height should be equal or greater than 100 meters (table 7).Since the basis of the hills have been recorded, it was sufficient to calculate the amplitude of the altitude in each of the hills, and classify the results leaving only those with greater value than 100 meters.
The function ZONAL STATISTICS [7.1] was used for calculating the amplitude of the altitude within the limits of the hills, obtaining values based on the DEM.Then, the hills with a height lower than 100 meters were  "dem_psaddle"; "range_hill" "third_01" Parameters "%DEM%" -"%h_saddle%"

Delimiting the upper third of the hills
Although a hill meets the average slope and the height criteria, the extent of the PPA comprises only the top third of it.So, the delimitation of PPAs considered the equivalent proportion of the higher third of the amplitude of each hill.
For each cell present within the area of the hills, the difference between the altitude in that cell (given by DEM) and the altitude of the saddle point of that hill (RASTER CALCULATOR [8.1]) were calculated (Table 8).In parallel, using the hill base as analysis zone, the amplitude of the altitude in each of these areas was calculated (ZONAL STATISTICS [8.2]).
The upper third was delimited by dividing the difference between the values of DEM and the altitude of the saddle point by the amplitude of altitude in each hill (RASTER CALCULATOR [8.3]).In the cells wherein this ratio was equal to or greater than 0.667, it was assigned a value "1" and the other with "0".Then, this raster was reclassified, leaving only the cells with value "1", which corresponded to the upper third of hills (RECLASSIFY [8.4]), illustrated in Figure 8. "ppa_rstr" Parameters ---Output "high_slope" "ppa_rstr" "ppa_hilltops" For this, the two raster files representing the hills with height greater than 100 meters and with slopes higher than 25 degrees (both containing only unit values) were multiplied (TIMES [9.1])(table9. A new raster was produced, whose cells identify the hills that meet the two criteria (height and average slope).This, in turn, was multiplied by the raster of upper thirds (TIMES [9.2]), resulting in Permanent Preservation Areas on the upper third of the hills.The result was then converted to the polygon vector format (RASTER TO POLYGON [9.3]), for ease of use in further analysis.

Delimitation of Permanent Preservation Areas
Once mapped, the hills that meet the criteria of height greater than 100 meters and average slope greater than 25 degrees, individually, all it took was to identify which of them meet both criteria.Thus, the upper third of them are classified as Permanent Preservation Area.

CONCLUSIONS
The methodological model presented in this work represents an opportunity in automatic mapping of PPAs, agreeing with the needs of the New Forest Code.This methodology could be used in the Rural Environmental Registration (CAR) system, which does not consider of PPAs on hilltops, to the present date.However, a field validation is needed to check the mapping results.This should assess both the quality of the methodological model as databases (DEMs) used.
Some issues must be defined by regulation in order to eliminate the subjectivity in mapping, namely the appropriate scale for mapping and the method for calculating the average slope.An offi cial database, validated by the government, should also be made available for this application.The resolution of the digital terrain model used could affect the results.A more specifi c approach to this matter will be made in future work.
The systematization of operations in Model Builder made possible the quick production of maps of PPAs on hilltops (Figure 9).The creation of large scale maps becomes feasible using the proposed method.
7), were converted to the polyline format (FEATURE TO LINE [3.1]).The new file, which contains the individual lines, was then used in the function ZONAL STATISTICS [3.2], where the lines were placed on the zone field and the DEM as input values.For each line, the maximum altitude of the DEM was then indexed, taking into account the region covered by each line separately.As result, a raster file where each rasterized line contained only a single value, its maximum altitude.Knowing the altitude where these lines intersect the saddle point, it was necessary to equalize this file to DEM by function RASTER CALCULATOR [3.3].The result was a new raster file, which contained only cells with the unit value, where the altitude of the DEM was equal to saddle point altitude, and zero for the others.Thereafter, the latest file was reclassified (Reclassify [3.4]), eliminating null values, leaving only the cells with value 1, which represent the saddle points.The result was then converted to integer type with function INT [3.5], and finally converted to the point format (RASTER TO POINT [3.6]).

FIGURE 1
FIGURE 1 Product generated in the first step: area boundary of each hill.

FIGURE 2
FIGURE 2 Product generated in the second step: the top point of the hills.

FIGURE 3
FIGURE 3 Product generated in the third step: saddle points.TABLE 3 Parameter settings in the functions used in the location of the saddle points.Names followed by asterisks represent only files where information was added in the attribute table.

FIGURE 4
FIGURE 4Product generated in the fourth step: limits of the hills with saddle points and top point data.

FIGURE 6
FIGURE 6Fifth step: hill base defined by the altitue of saddle point.

FIGURE 9
FIGURE 9 Example of a fi nal product.Map of PermanentPreservation Areas on hill tops.

TABLE 1
Parameter settings of the functions used in the delimitation of the boundary area of each hill.

TABLE 2
Parameter settings of the functions used to find the top points.Names followed by asterisks represent only files where information was added in the attribute table.

TABLE 5
Parameter settings in the functions used for delimitation of the hill base.

TABLE 4
Parameter settings in the functions used to determine the closest saddle points of the elevations.Names followed by asterisks represent only files where information was added in the attribute table.

TABLE 6
Parameter settings of the functions used to calculate the average slope of the hills.Names followed by asterisks represent only files where information was added in the attribute table.
FIGURE 7Product generated in the seventh step: hills that are higher than 100 meters.

TABLE 7
Parameter settings of the functions used to calculate the height of the hills.
FIGURE 8Product generated in the eighth step: upper thirds of hills.

TABLE 8
Parameter settings of functions used for defining the upper third.

TABLE 9
Parameter settings in the functions used for delimitation of Permanent Preservation Areas on hill tops.