Using Landsat Satellite Imagery to Determine Land Use / Land Cover Changes in Međimurje County , Croatia

This study evaluates Land Use and Land Cover (LU/LC) changes in Međimurje County, Croatia, from 1978-1992, 1992-2007 and 1978-2007 using Landsat satellite images. Spatial dynamics of LU/LC changes were quantified using three Landsat satellite images: MSS, TM and ETM+. Three series of the LU/LC maps of 1978, 1992 and 2007 were produced. Based on the Anderson classification system, LU/LC are classified as: water bodies, forest land, barren land, agricultural land and urban built-up land. Error matrices were used to assess classification accuracy. Significant buildings construction, construction of the reservoir lakes on the Drava River and construction of a highway were major drives of LU/ LC changes in Međimurje County over the study period. This study showed a continuous decrease in agricultural lands in Međimurje County. Consequently, the urban built-up area and the water area increased.


INTRODUCTION
It is well known that there are only few places on the Earth that are still in their natural state and that have not been affected by human activity in some way.These human activities result in significant land use changes at regional and local scales together with ecological, socio-economic and aesthetical impacts.Timely and precise information about change detection of Earth's surface is extremely important for understanding relationships and interactions between human and natural phenomena (Lu et al., 2004).Determining the effects of land use change on the Earth system, especially biodiversity, depends on the understanding of past land use practices, current land use patterns, and projections of future land use, as affected by human institutions, population size and distribution, economic development, technology, and other factors (Jingan et al., 2005).Viewing the Earth from space is crucial to the understanding of the influence of human activities and human impacts on land use changes over time.Information on LU/LC at regional scales derived from observations of the Earth from space provides objective information of human utilization of the landscape.Such information is required to support environmental policy, physical planning purposes and sustainable land use and land development.
Međimurje County, the northernmost and the most densely populated county in the Republic of Croatia has witnessed remarkable expansion, growth and development activities such as significant buildings construction, construction of two reservoir lakes on the Drava River and construction of a highway through the central part of the county.Such a rapid increase of land consumption and modifications on LU/LC changes resulted in lack of attempts to map and evaluate these changes.Therefore, the aim of this study was to produce LU/LC maps of Međimurje County and to identify and analyze general trends in LU/LC changes taking place in Međimurje County over a period of 29 years using Landsat Satellite Imagery and GIS based technique.

STUDY AREA
Međimurje County is a triangle-shaped county in the northern part of the Republic of Croatia and covers the plains between two rivers -the Mura and the Drava (Fig. 1).The elevation of Međimurje County ranges between 120 and 344 metres above sea level.The total area of Međimurje is 729 km 2 and it is the most densely populated county in Croatia.Almost 120,000 inhabitants live on 729 km 2 in 3 towns, 22 municipalities and cca.120 villages, amounting to the population density of 164.2 inhabitants/km 2 .The population of the Međimurje County accounts for 2.5 % of the total population of The Republic of Croatia.In the last 30 years there have been significant changes in environment such as: -Construction of two reservoir lakes on the Drava River -Lake Varaždin and Lake Dubrava -both built to serve the two hydroelectric power plants.-Construction of a highway through the central part of the Međimurje County.
-The rapid urbanization and industrialization, especially in the eastern part of the county.

METHODOLOGY, DATASETS AND SOFTWARE USED
The methodology used in the study includes the following stages: image pre-processing, the design of classification scheme, image classification, accuracy assessment and analysis of the LU/LC changes.The workflow diagram for the methodology used in this study is shown below (Fig. 2).
For the purpose of this study the multispectral Landsat satellite images of following were used (Tab.1): Landsat 4-MSS of August 23, 1978, Landsat 5-TM of August, 28 and Landsat ETM + of June 2007.The study area is contained within the Landsat path 204, row 028 for Landsat MSS and path 189, row 028 for Landsat TM and ETM + .All three series of the Landsat data images are created into global data set -Global Land Survey (GLS).The dates of the Landsat images were chosen to be as close as possible to the same vegetation season.The Landsat data images are georeferenced by the EROS Data Centre (Sioux Falls, South Dakota) to less than ½ pixel root mean square error, registered to Universal Transverse Mercator coordinates, zone 33N, WGS84 Datum.To ensure consistency among Landsat satellite images, all three series of the Landsat data images were re-sampled to a common nominal spatial grid of 30 meters resolution using nearest neighbour transformation to avoid altering the original values (Jensen, 1996;Yang and Lo, 2002).All three Landsat satellite images series are downloaded from the Global Land Cover Facility at the University of Maryland.
The following ancillary data (Tab.2) were procured: high resolution colour aerial photographs acquired between 2008 and 2009, black and white aerial photographs acquired in 1984 and 2002, topographic maps of 1:25,000 scales acquired in 1967 and 2001.
All procured ancillary datasets were used for "ground-truth" information required for classification and accuracy estimation of classified Landsat satellite images.The System for Automated Geoscientific Analyses (SAGA) GIS software (v.2.0.7) has been used for the preliminary data processing, extracting the study area, unsupervised and supervised classification and producing change detection maps.SAGA is a Free Open Source Software (FOSS).SAGA's analytical and operational capabilities cover geostatistics, terrain analysis, image processing, georeferencing and various tools for vector and raster data manipulation (Conrad, 2006).
Typical steps of pre-processing of the data are: geometric, radiometric and atmospheric corrections, re-sampling and sub-setting.Geometric correction of the data is critical for performing a LU/LC change detection analysis.Once the various data have been obtained, they have probably different projections, and the next necessary step is to project it.Thus, all obtained spatial data sets of a study belong to one single coordinate system.In this study the geometric correction of the Landsat satellite data has not been performed because the data had already been orthorectified and georeferenced to the UTM Zone 33N coordinate system.Accordingly, the Universal Transverse Mercator (UTM) projection, Zone 33N, Spheroid WGS84 and WGS84 Datum was used with that as the primary coordinate system.We have a different situation with ancillary data obtained from the State Geodetic Administration (SGA).Because the ancillary data are in three different coordinate systems to ensure consistency among the Landsat satellite images and ancillary data, we have to make geometric correction of the ancillary data.The geometric correction of all ancillary data that are used for this study has been done with SAGA and Proj 4. library.Proj 4. library work for raster and vector data and provide various projections for free definable cartographic parameters (Conrad, 2006).After the geometric correction of ancillary data, visual comparisons with a ground control point selected from topographic maps, as well as image-to-image comparisons showed that the images were geometrically correct.
In addition, the Landsat scenes and ancillary datasets are much larger than the study area and it is beneficial to reduce the size of the image files to include only the study area in the work process.This reduction of data is known as sub-setting, and this reduction eliminates the extraneous data in the file and speeds up processing due to the smaller amount of data to process.The sub-setting of all datasets has been done with SAGA as well.The primary classification scheme used for this study was based on the Anderson et al. (1976) land-use/cover classification system for a level one classification, along with the author's good knowledge of the terrain.The derived final LU/LC classification scheme utilized five classes: Urban/Built Up, Water, Barren, Agricultural and Forest (Tab.3).
Tab. 3 The classification scheme for the study Tab. 3. Klasifikacijska shema istraživanja

Class Descriptions
Urban/Built-Up Includes all residential, commercial, and industrial developments and transportation.
Water Includes all water bodies (rivers, lakes, gravels, streams, canals and reservoirs).

Agricultural
Includes all agricultural lands.

Forest
Includes all forest vegetation types (evergreen, deciduous and wetland), and non-forest vegetation futures that are not typical of forest (pasture grasslands and recreational grasses).
According to Meyer (1995), every parcel of land on the Earth's surface is unique in the cover it possesses.Land Use and Land Cover are distinct yet closely linked characteristics of the Earth's surface.Land Cover and Land Use are sometimes used as terms with the same meanings.Actually they have different meanings.The mixing of the concepts of LU/ LC has been present for at least the last 25 years (Anderson et al., 1976).The difference in LU/LC is acknowledged in many documents (Anderson et al., 1976;Campbell 1981, Di Gregorio andJansen, 2000;Cihlar 2000).Land cover is the physical material on the surface of the Earth.It is the material that we see in direct interaction with electromagnetic radiation and causes the level of reflected energy that we observe as the tone or the digital number at a location in an aerial photograph or satellite image.Land use, by contrast to land cover, is a description of how people use the land.(Fisher et al., 2005).Land cover categories include: grass, asphalt, cropland, water, snow, deciduous forests, bare soil, wetlands and pasture.Examples of land use include: agricultural land, urban and recreation areas, grazing and mining.
Training area selection is an important step in the classification.Training areas must be homogeneous and represent LU/LC classes.There are many ways to collect training data which include: collecting from field information (ground survey), on screen selecting training data as polygons, and on screen seeding of training data (Jensen, 2005).In this study, the obtained ancillary data, aerial photographs and topographic maps were used during the selection of training samples for "ground truth" information required for classification, and for accuracy assessment of achieved classified LU/LC maps.Also, the training areas were based on the author's good knowledge of the study area, and on visual comparison of the natural and false colour composite images.Areas depicting the various defined land covers were heads-up digitized from that aerial photography and topographic maps.At least 25 training areas for the defined classes (except for urban areas) were created for the 1978, 1992 and 2007 Landsat MSS, TM and ETM+ respectively, to train the images for supervised classification.It is obvious that collecting ground truth data through aerial photographs has an advantage over the conventional method of survey because of its speed, accuracy, and cost effectiveness.
Urban surface is heterogeneous and typically composed of a complex combination of features that are smaller than the spatial resolutions of the sensors (buildings, roads, grass, trees, soil, water), (Jensen, 2000).This characteristic of urban landscapes makes mixed pixels a common problem in medium spatial resolution data (between 10 to 100 m spatial resolutions) such as Landsat MSS, TM and ETM + .Such a mixture becomes especially prevalent in residential areas where buildings, trees, lawns, concrete, and asphalt can all occur within a pixel (Lu and Weng, 2005).For this reason, in this study, the approach recommended by Robinson and Nagel (1990), Harris and Ventura (1995), Reese et al., (2002) was used.With this approach, the urban areas, due to confusion with the bare soil, can be classified more accurately if it is done separately from the rural areas.Careful manual delineation on the topographic maps 1:25,000 and high area photographs was done around all urban areas of the Međimurje County.More than 150 polygons greater than 100 contiguous pixels, (with size of the surface larger than 9 ha) were obtained.Obtained urban areas were then clipped out of the image data and classified separately from other data with an unsupervised classification.Pixels classified as low or high density urban were masked out of the others Landsat data, while non-urban pixels were "put back" into the Landsat image data for further supervised classification.
In all classification methods the resulting images of the classification are probably speckled, which means that individual, isolated pixels of one class are surrounded by pixels of another class, and classified images usually have to be post-processed before integration into a geographic information system (GIS).For this reason, filters are usually applied to smooth or generalize the final classification images.With a majority filter, each pixel is recoded to the majority class of a neighbourhood defined by the filter (Stuckens et al., 2000).Majority filters reduce the "salt-and-pepper" effect typical for traditional per-pixel.In this study, a majority filter was used and applied for all three classification maps.An improvement of 3.47 % compared to a per-pixel classifier was achieved.
In a variety of studies, post-classification change detection has been considered to be one of the most suitable and commonly used methods for change detection (Singh, 1989;Mas, 1999, Jensen,1996;Peterson et al., 2004).Post-classification is a term describing the comparative analysis of spectral classifications for different dates produced independently (Singh, 1989).This method involves comparing two independent classified LU/LC maps from images of two different dates.In this study, the post classification technique was used based on a hybrid classification approach comprised of unsupervised classification based on "hill climbing" cluster algorithm and supervised classification based on the Maximum Likelihood Classifier (MLC) and introduction of human-knowledge from field experience to determine the amount of change over the study period.
The error matrix is the most widely used approach for image classification accuracy assessment and can be used to derive a series of descriptive and analytical statistics (Congalton and Mead;1983;Congalton, 1991;Congalton and Green, 1999;Smits et al., 1999;Congalton and Plourde, 2002;Foody, 2002;Liu et al., 2007).In this study, accuracy assessment was performed for the classified maps of all three time steps.Error matrices were used to assess classification accuracy using four measures of accuracy: overall accuracy, user's accuracy, producer's accuracy and Kappa statistics.The meaning of and calculation for these terms have been explained in detail in many studies (Congalton 1991;Smits et al., 1999;Plourde and Congalton, 2003;Foody 2002).The user's accuracy indicates the probability that a pixel on the image actually represents that class on the ground (Story and Congalton, 1986).It is calculated for each class by dividing the correctly classified pixels by the row total for that class.The producer's accuracy is defined as the probability of a pixel being correctly classified and is mainly used to determine how well an area can be classified (Story and Congalton, 1986).Over the past 15 years, Kappa statistics has become a standard part of evaluating classification accuracy.Kappa statistics is a discrete multivariate technique used to evaluate the accuracy of change detection and classification maps by measuring the agreement between the two images (Story and Congalton, 1986).Kappa statistics is a measure of the difference between the actual agreement between the reference data and an automated classifier and the chance agreement between the reference data and a random classifier.
Achieved results for the accuracy of all three years are summarized in Tabs.4, 5 and 6.For the accuracy assessment in this study, simple random sampling was adopted.A total of 146, 163 and 131 random pixels from the classified images of 1978, 1992 and 2007 respectively, without any consideration of informational class were selected.The pixels were verified using ancillary data.The overall accuracies for 1978, 1992 and 2007 were, respectively, 87.67 %, 88.96 % and 90.84 %.Kappa statistics were 80 %, 84 % and 83 %.User's and producer's accuracies of individual classes were relatively high, ranging from 74 % to 94 % which indicates a good agreement between thematic maps generated from images and the reference data.The achieved accuracies of classification turned out to be better than expected.

DISCUSSION
The most important LU/LC changes that have occurred in the Međimurje County over the study period from 1978 to 2007 were first analyzed by visual interpretation.Visual interpretation can give a general idea about the LU/LC changes.Using the Landsat satellite images natural and false colour, composite images were generated.Natural and false colour composites were chosen as a very easy, but very useful technique to extract visually information from Landsat satellite images and to achieve a general description of LU/LC changes.For Landsat TM and ETM + data, band 3 (red), band 2 (green), and band 1 (blue) were used to create a "normal colour" composite.Landsat MSS data do not have a blue band, and therefore they cannot be used to generate a true normal colour composite.However, the normal colour composite can be simulated using the following formula: Red=MSS2, Green=(MSS1*3 + MSS2) / 4 and Blue=MSS1.
For example, by comparing the three composite images from 1978, 1992 and 2007 we can see how the urban areas grow for the city of Prelog (Fig. 3).Further, we can see the expansion (Fig. 4) of the water surface due to the exploitation of the gravel over the years (images 1992 and 2007).On the time series of the three composite images (Fig. 5) we can see the construction of a reservoir lake on the Drava River built to serve the Čakovec hydroelectric power plant.All these changes detected by visual interpretation show us that there have been significant changes in LU/LC over the investigated time.
The LU/LC classification maps of 1978of , 1992of and 2007 were generated (Figs. 6, 7 and 8) were generated (Figs. 6, 7 and 8).Corresponding classification was also represented as a chart.The total area of each land cover class for the entire study was compared (Tab.7), using three classified images from 1978, 1992 and 2007.Međimurje County is predominantly a rural region.Tab.7 Land Use/Land Cover Distribution (1978, 1992 and 2007) Tab. 7. Distribucija načina upotrebe i pokrova zemljišta (1978., 1992. i 2007.)Land LU/LC changes maps for 1978-1992, 1992-2007and 1978-2007 were derived (Figs. 9, 10 and 11) were derived (Figs. 9, 10 and 11).The post classification technique was used to determine the amount of change over the study period.The total area of each land cover class for the entire study was compared (Tab.8) using three classified images from 1978, 1992 and 2007.Based on the results of this study, there is strong evidence that changes in the Land Use and Land Cover occurred in the Međimurje County during the last three decades.The changes in Land Use were determined for Međimurje County as follows: increase in water body areas by 10 % from 1978-1992 and by 4 % from 1992-2007; significant decrease in agricultural land by 13 % from 1978-1992, and very small increase from 1992-2007; increase in barren land by 107 % from 1978-1992, and decrease by 28 % from 1992-2007; constant increase in the built-up category by 10 % from 1978-1992, and by 4 % from 1992-2007; decrease in forest by 9 % from 1978-1992 but increase by 16 % from 199216 % from -200716 % from . Generally, from 197816 % from to 2007, urban area increased approximately 53 %, while agriculture decreased 12 %, whereas water, forest and barren land also increased by 14 %, 5 % and 49 % respectively.Achieved LU/LC change maps showed that urban land area increased 53 % due to urbanization that resulted from high-speed economic development, especially in the period from 1978 to 1992, water increased 14 % mainly due to the construction of reservoir lakes on the Drava River and land conversion in the gravels.Consequently, forest area decreased 9 % from 1978 to 1992, but increased from 1992 to 2007, especially along the border of the existing forest areas.Since the early 1990s, the share of agricultural population is constantly decreasing, and therefore there was more abandoned land which, after a while, becomes a forest.In total, the forest land in the study area has increased by slightly less than 1000 ha.Increase in urban and forest areas resulted in a substantial reduction in the area of agricultural land.Looking at the barren and agricultural land as a single class, total agricultural land decreased by 2 % during 1978-1992 and the significant decrease by 5 % was from 1992-2007.Total reduction in agricultural land over the study period was 7 %, indicating a permanent reduction of agricultural land in Međimurje County.

CONCLUSION AND RECOMMENDATION
The aim of this study was to produce LU/LC maps of Međimurje County and to identify and analyze general trends in LU/LC changes.An effective approach for making LU/LC maps and detecting changes in Međimurje County was presented without using a complex methodology such as fuzzy classification or object-based methodology.Basically, the use of the traditional per-pixel based classification in combination with a separate classification of the urban areas approach, a good knowledge of the research area, GIS tools, availability of quality ancillary data for the selection of training samples for "ground truth" information required for classification and for accuracy assessment of achieved classified LU/LC maps, demonstrate the potential of multi-temporal Landsat satellite imagery to provide LU/LC maps and to analyze changes in LU/ LC.Such results can be used as inputs to land management and policy decisions, to the development of sustainable urban land use planning decisions and for governmental or non-governmental organizations in a variety of purposes as up-to-date information.
The key findings of this study are as follows: Međimurje County is predominantly a rural region.Over half of the county is used for agriculture.In the study periods covered, the major identified Land Use and Land Cover classes of Međimurje County included agriculture, forest, barren land, built-up areas and water.
Significant buildings construction, construction of the reservoir lakes on the Drava River and construction of a highway were major drives of LU/LC changes in Međimurje County over the study period from 1978 to 2007.This study showed a continuous decrease in agricultural lands in Međimurje County.Consequently, the urban built-up area and the water area increased.
The remarkable changes in land use classes have occurred from 1978 to 1992, but the changes were not so significant during 1992 to 2007.
The remote sensing and classifications of the Landsat satellite imagery can be used as an economical and accurate way to produce Land Use and Land Cover change maps in the Međimurje County that can be used as inputs for regional analyses, formulating effective environmental policies and resource management decisions as useful up-to-date information.One of the most interesting questions when considering future work is the usefulness of high-resolution (spatial and spectral) satellite and airborne sensors, for detailed map and research of LU/LC changes.It would be very interesting in future work to use a combination of high resolution remote sensing imageries such as Quickbird and IKONOS, GIS and Landsat satellite imagery in order to improve the mapping, to detect changes in Međimurje County and for creating detailed Land Use and Land Cover maps.SAŽETAK Upotreba satelitskih snimaka Landsat za utvrđivanje promjena u načinu upotrebe i pokrovu zemljišta u Međimurskoj županiji u Hrvatskoj