Mapping of deforested areas in the city of Rio de Janeiro

Mapping of deforested areas in the city of Rio de Janeiro

(Parte 1 de 2)


Topic: Mapping of deforested areas in the city of Rio de Janeiro, approaching the changes from 1985 to 2014 using Landsat 5 and Landsat 8 satellite images

Prepared by:

Leandro Luiz Silva de França

December 2nd, 2014

ANNEX I – MatLab Script to Selection and Normalization of shaded vegetation18
ANNEX I – MatLab Script to Calculate the NDVI and determine vegetation areas21

Content ANNEX I – MatLab Script to determine deforested, unaltered and forested areas ................. 23

Satellite multispectral imagery has been widely used in geological and environmental studies. Despite of the low spatial resolution of a Lansat image (30 m), it is possible to extract vegetation information through NDVI (Normalized Difference Vegetation Index). In this work, a method of determining deforested areas based upon the threshold the NDVI to specify vegetation is presented.

The objective of this work is to map the deforested areas, as well the areas where the vegetation has remained unaltered and forested areas in the city of Rio de Janeiro considering the time’s interval from 1985 to 2014. The results of the recognition of each specific area are shown in a map.

Multispectral imagery has become a good beginning in the study of remote sensing. Even though the low spatial resolution interferes the accuracy of the determination of the exact identification of features in the terrain, those images can be useful to determine great patterns depend on the scale of the study area. That means for vegetation measurement those imagery may be used to analyze great modifications.

The NDVI is one of the most successful and fast way to detect live green plant canopies, because they scatter more near-infrared radiation than others features. The calculus of the NDVI for each individual measurement is the ratio between the difference of near-infrared (NIR) and red (R) measurements by the sum of NIR and R measurements. As shown in the formula below:

The NDVI values can vary between -1.0 and 1.0, and for vegetation coverage they tend to be greater than zero.

However, hill’s shaded areas can interfere in the NDVI’s calculus (see Figures 1).

The Figure 2 presents a NDVI image with pixels values greater than zero in purple. We expect that all vegetation area has NDVI greater than zero, however this does not happen. The shaded vegetation has NDVI less than zero.

Figure 1: Landsat RGB image of a rough relief area

Figure 2: Overlay between RGB image and NDVI greater than zero in purple. Shaded vegetation has NDVI less than zero

In the case of rough relief it is possible to estimate the pattern of shaded areas.

Based on this, it was used a process to identify and to eliminate the shaded aspect normalizing the pixel’s values in order to match to a non-shaded area with similar pattern.

The mentioned operation consists of previously take samples of shaded vegetation and samples of non-shaded vegetation (as shown in Figure 3). After that, the histogram is calculated (Figure 4 and Figure 5).

Figure 3: Sample areas for vegetation (green) and shaded vegetation (red)

Figure 4: Histograms for the shaded areas describe low values for each band

Figure 5: Histograms for the non-shaded areas show greater pixels’ value for all bands. We can note the significant values for NIR band (as expected)

The normalization enforces the histogram of the shaded vegetation be similar the histogram of the non-shaded vegetation (Figure 6). The calculus can be performed for each band (only in the shaded areas using as parameter the non-shaded vegetation distribution) through the following formula for each pixel’s value of a band (B),

() .( )shaded new vegetation vegetation shaded where, shadedµ and shadedσcorrespond the mean and standard deviation of shaded area.

And 3.vegetationvegetationµσ+ and vegetationσ correspond the mean and standard deviation of non-shaded area.

The Figure 6 shows the normalized histograms and the Figure 7 is the output after the image processing. The pixels with the same pattern of the shaded vegetation areas (mathematically when the pixel’s value is between the interval [

3.vegetationvegetationµσ−, 3.vegetationvegetationµσ+] , that is 9% of a Gaussian distribution) was transformed by normalization.

Figure 6: Histogram of shaded areas after the normalization process

Figure 7: Most of the vegetation with shaded areas was transformed in order to have similar pattern of non-shaded vegetation.

The Figure 8 shows the NDVI image of the non-normalized image. It is evident that the shaded areas interfere in the calculus of the NDVI. The second image (Figure 9) corresponds to the NDVI image of the normalized image. So, we can note that the previous problem was solved.

Figure 8: NDVI image of the non-normalized image

Figure 9: NDVI image of the normalized image

Nevertheless, we clearly observe (Figure 10) that the normalization procedure causes some side effects. In other words, some pixels in the region of water have the same pattern as the pattern calculated for shaded areas.

Figure 10: Side effects because of the normalization procedure (some pixels in the water got the same pattern of the vegetation

In order to fix the error of the normalization analyzing the four bands, an additional image called slope should be used. The slope (Figure 1) is an image which is generated from the digital elevation model (DEM) in which each pixel corresponds to the angle between the Earth surface and the flat horizon. This property should be applied as criteria to distinguish the pixels in the water of the pixels in the shaded areas.

Figure 1: Slope image and shaded areas (red)

The Figure 12 depicts the slope distribution for the shaded areas. As awaited, most of the values are in the interval of 20° and 50°.

Figure 12: Slope distribution

Therefore, in process of selection of the pixels with shaded vegetation areas for normalization, the slope patter is also verified in order to correctly match the sample pattern with all pixels with similar characteristics obeying this new criteria combined with the old.

Figure 13: Result of the normalization process using the slope as additional criteria to select the pixels (no more green pixels in the water)

The deforested, forested and unaltered areas of vegetation can be determined through Boolean overlay. The Figure 13 is an example that explains the logical operations to classify each case.

Figure 14: Boolean overlay Mathematically,


(_)(_)(_)deforestedoldvegetationnewvegetationnewvegetation=∪− (4) (_)(_)(_)forestedoldvegetationnewvegetationoldvegetation=∪− (5)

This section is dedicated to describe each step performed to achieve the research’s objective. The steps are:

1. Determine the study area and download its shapefile. In this case the shapefile of Rio de Janeiro's city boundaries, provided by IBGE (Instituto Brasileiro de Geografia e Estatística – Brazilian Institute of Geography and Statistics). 2. Download two satellite images of the study area in USGS Global Visualization

Viewer. The first is a Landsat 5 image which acquisition date is August 05, 1985. And second is a Landsat 8 image which acquisition date is January 25, 2014. 3. Select the bands to be worked. It was selected the blue, green, red and infrared band of each image. Respectively, for Landsat 5: B1, B2, B3 and B4, and for Landsat 8: B2, B3, B4 and B5. 4. Download the DEM for the study area. A DEM with 90 m of spatial resolution is provided by EMBRAPA (Empresa Brasileira de Pesquisa Agropecuária – Brazilian Company of Agriculture Researches). 5. Generate a slope image from the DEM. 6. Clip all images to same extent of the shapefile of Rio de Janeiro’s Boundaries and Resample then taking the Landsat 8 image as reference. 7. Create polygons for samples of shade and vegetation areas for 2014 and 1985 images. 8. Use the algorithm (Annex I) to normalize the areas with shaded vegetation for 2014 and 1985 images. 9. Calculate the NDVI and the threshold value which determines vegetation for 2014 and 1985 images (Annex I). For both images, the threshold is determined accordingly to samples in the same location and same green pattern in the two years, that is, grass and small trees. 10. Use map algebra (with Boolean overlay) to detect the unaltered vegetation coverage, deforest area and forested area (Annex I).

The calculated thresholds correspond to the minimum NDVI values to consider as vegetation and they were obtained by heuristic. That means some areas was arbitrary chosen considering that they virtually have not changed (yellow polygon in Figure 15 and Figure 16). Thereafter, it was calculated the NDVI for all pixels inside those polygons for each image. From those values, the mean and standard deviation was determined.

The adopted threshold values in this research are presented in the Table 1. The threshold is calculated by the equation (6). The condition above causes that more than 9% of the pixels inside the sample polygon are considered vegetation.


Tab 1: Threshold values

Image Mean Std. Deviation Threshold 1985 0.4413 0.0829 0.1925 2014 0.4122 0.0235 0.3418

The greater variation for the 1985 image is probably related to inferior quality of the sensor (Landsat 5) which counted on 8 bits of radiometric resolution, differently of the new sensor (Landsat 8) which counted on 16 bits of radiometric resolution.

Another likely reason for the greater variation is the different acquisition period of the year. The 2014 was obtained in January (rainy and hot period); on the other hand the 1985 was obtained in August (dry and cool period). So, the weather conditions can interfere the NDVI because they interfere the plants’ health condition and hence the amount of chlorophyll and cell structure of the leaves (responsible to absorb the light and reflect the infrared radiation).

Figure 15: Vegetation sample (yellow) to be consider the threshold pattern (1985)

Figure 16: Vegetation sample (yellow) to be consider the threshold pattern (2014) Figure 17: Vegetation change (1985 – 2014)

The Figure 18 is a full extent view of the study area, where the red places correspond to deforested areas, the blue areas are the unaltered vegetation and the green areas are the forested.

In a short analysis, we can note that the west part of the city suffered more deforestation than the east part. On the other hand forested areas were more presented in the north and east of the city.

Other evident and worrying fact is that the deforestation areas are mainly concentrated around the unaltered vegetation. It can mean that the deforestation is going to be increasing in the direction of the preserved areas.

Figure 18: Unaltered area (blue), deforested area (red) and forested area (green)

The objective of this work is materialized in the vegetation change map which shows the product of the methods and procedures used in this research.

Although the limitation of the spatial resolution and sometimes the quality of the images, it is possible to obtain good results in order to accomplish environmental analysis. Of course that it will depend upon the scale of the study area.

A suggestion for future researches is to test all presented procedure in other images, studying other variables, for example weather conditions. So, the threshold values should be examined further.


USGS (Landsat 5 History):

INPE (Divisão de Geração de Imagens):

PROCESSAMENTO DIGITAL (Landsat-8: Novas Combinações de Bandas e Informações Técnicas):


ANNEX I – MatLab Script to Selection and Normalization of shaded vegetation

% Experimento I - Normalizando pelas 4 bandas da imagem + Slope

% parametros (largura do intervalo e 9% - 3 desvios-padrão) largB = 3; largG = 3; largR = 3; largIR = 3; larg_slope = 3;

% 1 - Abrir imagem RGB e verificar áreas de amostras

% Abrir imagem GEOTIFF e separar em bandas filename = uigetfile('*.tif','Selecione um Arquivo GeoTiff Banda B'); [img, refmat, bbox] = geotiffread(filename); B = img(:,:,1); G = img(:,:,2); R = img(:,:,3); IR = img(:,:,4); filename2 = uigetfile('*.tif','Selecione um Arquivo GeoTiff Slope'); [slope, refmat2, bbox2] = geotiffread(filename2); origem = [bbox(1,1) bbox(2,2)]; % Coordenadas da origem da imagem escala = abs(refmat(2,1)); % Resolução da imagem tamanho = size(img); clear img bbox refmat bbox2 refmat2

% 2 - Abrir Shapefile e Preencher a estrutura com os polígonos das amostras Shape = shaperead(uigetfile('*.shp','Selecione um Arquivo Shapefile de polígonos para Amostra de SOMBRA')); num_poly = size(Shape,1); for k=1:num_poly tam = size(Shape(k).X,2); Xpoly = zeros(tam-1,1); Ypoly = zeros(tam-1,1); for j = 1: tam-1

% Translação com origem no primeiro pixel/ mudança de escala (unidades em px) Xpoly(j,1) = (Shape(k).X(1,j)-origem(1))/escala; % Translação com origem no primeiro pixel / inversão do eixo Y / mudança de escala (unidades em px)

Ypoly(j,1) = -1*(Shape(k).Y(1,j)-origem(2))/escala; end s(k).Polygon = [Ypoly Xpoly]; % Coordenadas de imagem end clear Shape Xpoly Ypoly k j

% 3 - Para cada polígono verificar os pixel que estão dentro de sua área

AMOSTRA = false(size(R)); for k = 1: num_poly

Poly = s(k).Polygon; % Determinando o retângulo que contem o polígono Xmin = floor(min(Poly(:,2))); Ymin = floor(min(Poly(:,1))); Xmax = floor(max(Poly(:,2))); Ymax = floor(max(Poly(:,1))); [X,Y] = meshgrid(Xmin:Xmax, Ymin:Ymax); IN = inpolygon(X-0.5,Y-0.5,Poly(:,2),Poly(:,1)); % -0.5 (centro do pixel dentro do polígono)

AMOSTRA(Ymin:Ymax,Xmin:Xmax)=IN; end clear Poly s Xmin Ymin Xmax Ymax X Y IN

vectorB(1,cont) = B(k);vectorG(1,cont) = G(k);
vectorR(1,cont) = R(k);vectorIR(1,cont) = IR(k);

% 4 - Gerar vetor de amostra com o pixels selecionados tam = size(AMOSTRA,1)*size(AMOSTRA,2); tam_vector = sum(sum(AMOSTRA)); % Alocando memória para os vetores vectorB = zeros(1, tam_vector); vectorG = zeros(1, tam_vector); vectorR = zeros(1, tam_vector); vectorIR = zeros(1, tam_vector); vectorSlope = zeros(1, tam_vector); cont = 1; for k = 1:tam if AMOSTRA(k) vectorSlope(1,cont) = slope(k); cont = cont + 1; end end

% 5 - Calcular a média e desvio-padrão para cada banda EstatBsombra = [mean(vectorB) std(vectorB)]; EstatGsombra = [mean(vectorG) std(vectorG)]; EstatRsombra = [mean(vectorR) std(vectorR)]; EstatIRsombra = [mean(vectorIR) std(vectorIR)]; EstatSlopeSombra = [mean(vectorSlope) std(vectorSlope)];

% 6 - Varrer cada banda determinando quais pixels estão dentro do intervalo % de aceitação para ser considerados similar ao da amostra.

IntervaloSlope = [ EstatSlopeSombra(1)-larg_slope*EstatSlopeSombra(2),

IntervaloB = [ EstatBsombra(1)-largB*EstatBsombra(2), EstatBsombra(1)+largB*EstatBsombra(2)]; IntervaloG = [ EstatGsombra(1)-largG*EstatGsombra(2), EstatGsombra(1)+largG*EstatGsombra(2)]; IntervaloR = [ EstatRsombra(1)-largR*EstatRsombra(2), EstatRsombra(1)+largR*EstatRsombra(2)]; IntervaloIR = [ EstatIRsombra(1)-largIR*EstatIRsombra(2), EstatIRsombra(1)+largIR*EstatIRsombra(2)]; EstatSlopeSombra(1)+larg_slope*EstatSlopeSombra(2)];

(Parte 1 de 2)