DIGITAL AERIAL IMAGES LAND COVER CLASSIFICATION BASED ON VEGETATION INDICES

Knowledge of how land cover has changed over time improve assessments of the changes in the future. Wide availability of remote sensed data and relatively low cost of their acquisition make them very attractive data source for Geographical Information Systems (GIS). The main goal of this paper is to prepare, run and evaluate image classification using a block of raw aerial images obtained from Digital Mapping Camera (DMC). Classification was preceded by preparation of raw images. It contained geometric and radiometric correction of every image in block. Initial images processing lead to compensate their brightness differences. It was obtained by calculating two vegetation indices: Normalized Difference Vegetation Index (NDVI) and Green Normalized Vegetation Index (gNDVI). These vegetation indices were the foundation of image classification. PCI Geomatics Geomatica 10.2 and Microimages TNT Mips software platforms were used for this purpose. Key worDS: image classification, land cover, digital aerial image, vegetation index, remote sensing. Piotr Dzieszko, Institute of Geoecology and Geoinformation, Poznań, Adam Mickiewicz University, ul. Dzięgielowa 27, 61-680 Poznań, Poland macieJ DzieSzKo, Piotr DzieSzKo, Sławomir Królewicz, Jerzy cierniewSKi


Introduction
Land is the main resource controlling primary productivity in terrestrial ecosystems (Darwin et al. 1996).Changes in land cover and land use af-fect the global systems (e.g., atmosphere, climate and sea level) (Meyer & Turner 1992).Information describing current land cover is an important input for planning and modelling, but the quality of such data defines the reliability of the simulation outputs (Townshend 1992;Belward 1996).The relationship between land use and land cover are not always obvious.Land cover classes can support multiple uses.Land use can be defined as the human use of land while land cover can be defined as the biophysical state of the earth's surface (Turner et al. 1995).Knowl-Knowledge of how land cover has changed over time improve assessments of the changes in the future.Easy availability of remote sensed data and relatively low cost of their acquisition make them very attractive data source for Geographical Information Systems (GIS).Aerial and satellite images present reliable data which can be easily interpreted by an analyst.The analysis of remote sensed data takes relatively more time than data acquisition.That is why in last years it can be noticed finding a solution which will give an opportunity to interpret remote sensed data in more and more automatic way.A good example of such automation can be fast updating of topographic databases, creating digital elevation models and land cover maps.With recent advent of very high resolution satellite imagery (HRSI), such as IKONOS, QuickBird, and OrbView images, great efforts have been made in the application of these remote sensing images in environmental studies.HRSIs have been applied widely in ur-HRSIs have been applied widely in ur-HRSIs have been applied widely in urban land-cover mapping (Im et al. 2008;Thomas et al. 2003;Wulder et al. 2008;Lu & Weng 2009).But when we want to map land cover characteristics in smaller scale, aerial images still have and would have a significant value.
It is impossible to create land cover map using available aerial images without previous analysis of them.But interpretation of images block is a complex case.Every image in block has different exposure and observation conditions.Terrain shape and land cover types have also significant influence on every image.In many cameras there is also vignetting effect which causes changes in exposure intensity in the image surface.These are the reasons why geometric and radiometric correction is obligatory in images analysis (Konecny 2003).It leads to eliminate potential image errors and to image normalization.
The main goal of this paper is to prepare, run and evaluate aerial image classification using a block of raw images obtained from Digital Mapping Camera (DMC).Classification was pre-ceded by preparation of raw images.It contained geometric and radiometric correction of every image in block.Initial images processing lead to compensate their brightness.It was obtained by calculating two vegetation indices: Normalized Difference Vegetation Index (NDVI) and Green Normalized Vegetation Index (gNDVI).These vegetation indices were the foundation of image classification.PCI Geomatics Geomatica 10.2 and Microimages TNT Mips software platforms were used for this purpose.

Study area
The study area is Świętoszów military training ground.It is located in Dolnośląskie Province.It is the biggest military training ground in Poland.It has about 38,000 hectares.Military training ground is traversed by many tank roads.There is a lot of military objects within its boundaries.Figure 1 shows the location of the study area.
The study area is located in the extent with higher precipitation values than observed in surrounding areas (about 600 mm annual) which makes that forest species composition more similar to upland forests.It is mainly beech forest with fir and spruce.

Input data
The block of digital aerial images were used to achieve the goal of the work.131 images collected in five strings were available for this elaboration.12 images were selected to analysis and they cover the centre of aerial raid area.The set of selected images contained 38, 39, 40 and 41 from third string, 63, 64, 65, 66 from fourth string and 81, 82, 83, 84 from fifth string.Strings have direction from the east to the west.The spatial resolution of images was 0.4 m.Data set contained images in five spectral channels: panchromatic, red, green, blue, and infrared.Digital processing of images was carried out with PCI Geomatics Geomatica 10.2 and Microimages TNT Mips software platform.
DMC was used to collect images.DMC was constructed in 2001 by Z/I Imaging.DMC is multi-module camera which allows to obtain 5 cm spatial resolution.It is possible thanks to using high quality CCD matrix with very high spectral resolution (12 bits for every channel).Focal length of the camera equals 120 mm and the available matrix size is 7,680×13,824 pixels.Images were collected from height of about 4,100 meters (Intergraph 2010).
In image classification, the use of ancillary data is important and effective especially in the preclassification, stratification and classifier modification (Harris & Ventura, 1995;Williams 2001).What is important, for digital remote sensing, ancillary data must be incorporated into the analysis in a structured, formalized manner (Campbell 2007).Unfortunately, for study area no ancillary data was available.

Methods
The work for obtaining the goal which was the image land cover classification was established according to workflow presented in Figure 2.

Digital Elevation Model (DEM) creation
To orthorectify aerial images high resolution and good quality DEM is necessary and crucial.Presented in Figure 3 10 m spatial resolution DEM was delivered with images but we came to conclusion that such a spatial resolution cannot guarantee enough precision for image analysis.Images spatial resolution was 0.4 m so spatial resolution of DEM was, in this case, more than 20 times lower.So another DEM was calculated with spatial resolution 1.6 m using photogrammetric methods on the basis of aerial images.To create DEM, set of 30 images was used.It con-  tained 12 images selected for the analysis and images neighbouring with them to obtain fine results on the selected images' edges.
Exterior orientation data was collected from Global Positioning System / Inertial Navigating System (GPS/INS) report stored in the aerial raid report.Firstly, DEMs for epipolar images were calculated.Next, DEMs were mosaicked in strings and after that, in whole area.Such an attitude allowed to avoid errors on the contact area of images.After all, final DEM was filtered with median filter in Microimages TNT Mips.Results of filtering are showed in Figure 4. Final DEM used in the analysis is presented in Figure 5.

Images orthorectification
When DEM for study area had been created images orthorectification was carried out.Green, red and infrared spectral channels images were orthorectified because exactly these images were necessary to calculate NDVI and gNDVI indices.PCI Geomatics Geomatica 10.2 was used for orthorectification.Figure 6 presents one image before and one after orthorectification.Running orthorectification with using digital elevation model with 1.6 m spatial resolution allowed for obtaining images more appropriate for vegetation indices calculation than images orthorectified using 10 m DEM.Their precision, especially in forested areas was much higher.Figure 7 presents differences between the same area orthorectified with 10 m and 1.6 m DEM.Using higher resolution DEM forested areas are better mapped than using lower resolution DEM.

Vegetation indices calculation
Images land cover classification was carried out on the basis of two vegetation indices (NDVI and gNDVI).These indices describe vegetation development state, its health state, they also describe biomass amount and participation of soil and vegetation in spectral reflectance from the unit area.NDVI and gNDVI indices were calculated in Microimages TNT Mips software using SML language script for automation the calculating process according to equations: NDVI (Normalized Difference Vegetation Index) , (1) gNDVI (Green Normalized Difference Vegetation Index) , (2) where: NIR -data from infrared spectra channel, R -data from red spectra channel, G -data from green spectra channel.
SML is a general-purpose, modular, functional programming language with compile-time type checking and type inference.SML script used for vegetation indices calculating is attached as an Appendix 1.The main task of the this script was to carry out the algebra operations according to equations which describe NDVI and gNDVI indices.Script used data from appropriate spectral channels and generate new images with vegetation indices values.In such a way 24 new images were obtained, 12 images for every index.Normally NDVI and gNDVI indices values spread from -1 to 1. Their ranges in this work were set from -10,000 to 10,000 to avoid necessity of decimal point using.
NDVI and gNDVI images had to be mosaicked to create a homogeneous surface with normalized vegetation values.Indices' values in single images were dependent on grid brightness in each spectral channel.Mosiacking allowed to solve this problem.
Mosaicking of images was carried out in Microimages TNT Mips.On overlapped area indices values were calculated with feathering function.This function uses for image normalizing weighted average and weight value is higher when grid is further from the image's edge.It gives much more better results than using simple arithmetic mean.The difference between results of these two methods presents Figure 8.
Figure 9 shows final mosaics for both indices: NDVI and gNDVI calculated for the study area.

Supervised classification of land cover
The very first step of the image classification was initial land cover class definition.Using aerial images for study area, general overview of the area was performed.The main four land cover classes were in this step identified.Set of initial classes contained (A) sands without vegetation cover, (B) grasses and scrub class, (C) heather class and (D) forested areas.Occurrence of these land cover types in study area is presented in Fig- ure 10.Table 1 shows vegetation indices values for first identified classes.Indices' values for land cover classes were calculated on the basis of correlation graph between NDVI and gNDVI mosaics.Correlation graph with simultaneous view in  Microimages TNT Mips allowed for land cover classes definition.Correlation graph is presented in Figure 11.After initial four classes definition further grid location on the correlation graph was carried out.Grid location on the graph was basis in process of unclassified grids assigning into destination classes.During this process it turned out that delimitation of another classes is necessary.Whole analysis lead to 12 classes definition.They were described by 25 intervals for every vegetation index which equals 50 intervals taken into account in the analysis.Determined classes with their intervals values are contained in Table 2.In the Figure 12, test areas with determined classes are presented.

Land cover classification performance
After precise interval vegetation index values determining for each land cover classes, classification could be performed.SML script was used for classification automation.The SML script code is contained in Appendix 2. The main goal of this script was to retrieve values of both vegetation indices, create a new raster file and assign to it values according to developed rules contained in Table 2.These rules were described by minimum and maximum values of intervals for both indices.For example, grid with NDVI value 4,000 and gNDVI value 2500 was assigned for coniferous forest class marked number 8. In the final classification map this grid obtained value 8 according to class where it was classified.Iteratively, for a every grid in such a way, grids were classified.In special cases new assigned grid value equaled 0. It meant 'no classification'.On the basis of developed script 85% of grids were classified.Remaining 15% of grids could not be classified to any class because development of full classification rules for such a big raster (about 560 billion of grids) was not possible in available software.What is also worth noticing, orthorectified images are north directed.Pixels located within the image frame were also assigned as 0 value.It explains such a big participation of not classified grids in final map.

Results
Figure 13 shows final classification map created on the basis of vegetation indices.New colour ramp was created for this map for better presentation of classes.New colour ramp with classes description is contained in Table 3.When classification was carried out it was necessary to evaluate it.It was done using 5 test areas and identifying occurring problem within them.This method was based only on remotely sensed data because field verification was impossible for study area.Civilians cannot enter area of interest chosen for this work.So classification precision was evaluated only using digital images interpretation.Especially useful was in this case orthophotomap delivered by courtesy of Forest Management Office.
Three main classification problems were identified and investigated: Shadows -they occurred often in final classification raster map.Shadows cannot be classified as a separate class but they were distinguished for better classification of remaining classes.Shadows mainly occurred on forested areas.That is why on the generalization step they should be assigned to forest class.This problem is presented on the Figure 14.
Not classified areas -these areas were mainly observed in the image frame area.In the main classification raster map content they occurred because vegetation indices values intervals were not closed in some cases.The effect of this was that some grids could not be assigned to any class.Problem of occurring not classified grids presents Figure 15.
Misclassification -which means assigning some grids to a wrong class.In every situation it was caused by very similar spectral characteristics of both land cover types.It was then impossible to distinguish such a regions on correlation graph.Such a cases were observed in following examples: Classification of bright sands and tarmac.Very interesting example concerns the old surface of A18 highway.New lane of A18 highway was assigned correct to tarmac class and old lane was assigned wrongly as a bright sand.It is shown on Figure 16.
Distinction of mixed and coniferous forest.What is worth stressing, using this method classification of forest class was good and did not cause any problems but it was sometimes difficult to identify forest type.Mixed forest normally   18.
Problem with assigning vegetation class on wetlands.In some areas wetlands with veg-  etation were classified as water.This problem presents Figure 19.
Every method has its advantages and disadvantages.These identified problems did not decrease quality of classification in authors' opin-ion.In the process of intervals filling, 12 classes were distinguished.Using developed method in future, it is good to decrease number of classes which can cause a solution for some identified classification problems.The main goal of this work was achieved and land cover classification was carried out and finished with success.The efficiency of classification was evaluated for 85% but the majority of misclassified grids are contained in image frame.
Calculating vegetation indices lead also to normalize images block which is definitely advantage of this method because in the step of image preparation none of other normalization technique is required.

Discussion and conclusions
In this developed method of digital image land cover classification some things need to be stressed.The work revealed that digital image preparation is time-consuming and complicated process.Special attention need to be put for the process of digital image normalization and one has to be aware to what consequences it may lead.
First step of the work revealed that precise DEM has a significant influence on whole analysis.It increases the quality of orthorectification and what is more, higher orthorectification quality allows to better calculation of vegetation indices.It is important in this case because vegetation indices were the basis of whole classification analysis.
Two used vegetation indices (NDVI and gNDVI) in complex way describe study area.The advantages of these indices were proved in this paper, especially their influence on image normalization and their efficiency in image classification.
The use of feathering method during images mosaicking revealed its advantages and allowed to skip the process of linear trend removing.
Classification itself was a little time-consuming mostly because intervals for vegetation indices for specific classes had to be identified as a rectangles.Analysing vegetation indices correlation graph it appeared that cloud point of specific class could be identified as an ellipse.It could better imitate classes relationships but using ellipses would increase time of analysis.That is why authors resigned from using ellipse's shape in the analysis.
Analysis should be more precise if there would be possibility of classification verification in the field.Because of the military character of study area, such a verification could not be performed.So land cover types were determined on the basis of image interpretation.It has influence on classification quality because some grids were assigned to incorrect class.The decrease of class number should increase classification quality.
The main advantage of this method is the possibility of land cover map creating with relatively short time and without time-consuming field works.

Fig. 3 .
Fig. 3. 10 m spatial resolution Digital Elevation Model delivered by Forest Management Office with input data.

Fig. 4 .
Fig. 4. Comparison of Digital Elevation Models.Before filtering by median filter (A) and after filtering (B).

Fig. 5 .
Fig. 5. Digital Elevation Model with spatial resolution 1.6 m created for the study area.

Fig. 7 .
Fig. 7. Differences between the same area, one orthorectified with 10 m DEM (A) and another orthorectified with 1.6 m DEM (B).

Fig. 10 .
Fig. 10.Four main land cover classes in study area.

Fig. 13 .
Fig. 13.Final classification raster map on the basis of NDVI and gNDVI indexes.

Fig. 14 .
Fig. 14.Shadows on the classification raster map.They are marked with white colour on the final raster map.

Table 1 .
Summary of vegetation indexes values for initial four land cover classes.

Table 2 .
Land cover classes determined on the basis of vegetation indexes for study area.

Table 3 .
Color ramp and classes description for final classification raster map on the basis of NDVI and gNDVI indexes.
oped for forest type distinction did not check out in 100%.This problem is shown on Figure17.

Misclassification of dense heather and conif- erous forest.
Dense heather has very high values of NDVI and these values overlapped with NDVI values for coniferous forest.It causes problem with misclassification of these two classes.This problem is presented on Figure