1. Introduction
The use of remote sensing data for archaeological research is not a recent development and can refer to any method that detects anthropogenic features or artifacts without direct contact with that material (Parcak 2009). This definition encompasses two different kinds of research, ground-based methods, such as magnetometry, and non-land-based methods, such as aerial photography. Ground-based remote sensing relies on things like ground penetrating RADAR or changes in magnetic field, as detected from background, to detect artifacts (206, Kvamme 2006). These techniques give a detailed example of what is buried in the ground, mapping out subsurface artifact scatters or buried landscape features without excavation (32, Slater et al. 1999), and in enough detail to map out individual hearths (Marshall 1998).
Non-land-based, usually aerial or satellite, remote sensing explores archaeological material with a “bird‘s eye view” (Parcak, 2009). However, this method neglects most artifacts, such as potsherds, since they are small in comparison with the resolution of the images. Therefore, archaeological uses of remote sensing focus on the detection of prehistoric built environments, such as tells or barrows, rather than individual artifacts (Parcak 2009, Saturno et al. 2007, 124 Bevan and Conolly 2004). The use of this data beyond the detection of possible sites, associated with landscape features, represents the limit of archaeological uses of remote sensing. Outside of archaeology, other disciplines use satellite remote to track changes in plant growth or to monitor of wildfires (Lasaponara 2006, Kuemmerle et al. 2006).
This paper explains a new use of remote sensing in archaeology, one that follows the footsteps of ecology. The goal of this new method is to use free or low-cost satellite remote sensing data in order to plan archaeological survey routes. While the application of this may seem redundant for current archaeological remote sensing projects (Saturno 2007), many archaeological projects utilizing remote sensing have focused on areas with little or no canopy (64, Parcak 2009) allowing many studies to detect archaeological features with the naked eye (Saturno 2006). However, many archaeological studies struggle with visibility and accessibility because high-density vegetation covers up artifacts and can limit survey in an area (16, Given 2004; 174, Llobera et al 2010). The use of remote sensing data to plan survey routes before entering the field would allow archaeological surveys to focus on areas that are easy to survey and maximize in field efficiency. For the purposes of this project, “easy to survey areas” are defined as places that would have little or no vegetation. While many projects ground-proof their findings to assess their validity (2281, Calvao and Palmeirim 2004), one of the aims of this project is to see the potential of this method as a planning tool. Therefore, ground-proofing before implementing the method would defeat the purpose.
This paper uses Normalized Difference Vegetation Indices, explained in detail below, in order to locate areas that have photosynthetic patterns or low NDVI and evaluates the usefulness of this method. Areas with low NDVI are assumed to have low vegetation density and would be easy to walk through. Utilizing the MODIS MOD13Q1 product, this study explores seasonal changes on Corsica for the last five years comparing changes in NDVI from winter to summer and from summer to winter. Areas that have seasonal changes in photosynthetic level, from highly photosynthetic to scantily photo synthetic or vice versa, represent places where the majority of plants have an “off season” and are less active. It is during this “off season” when survey could be executed and combined with survey of other low vegetation density regions.
1.2 Background on Study Area
The Mediterranean is an area of intense archaeological study (Given, 2004) and in some areas, ground-walking survey records so much material that clicker counters are the easiest way to count artifacts consistently (46, Fentress 2000). However, these high-density sites do not rule the Mediterranean. While a rather populous and large island, archaeological study of Corsica has not lead to artifact rich sites and therefore the island has not had the intense study of other regions (169-170, Llobera et al., 2010). Renewed study has begun with the hopes of systematic survey of the island. These studies, executed over multiple years, explored seven communes on Corsica, though survey for each commune did not exceed 10% of the area (Table 3. 175, Llobera et al, 2010). Before archaeological survey occurred municipal permits were acquired for each commune but the real right to survey rested with the landowner, or owners, rather than with government permits (173, Llobera et al, 2010). These logistical inconveniences did not compare to the problems posed by the flora of the island. The French word maquis describes scrublands found all over the Mediterranean. The term does not refer to a single plant, but instead to a mixture of plants, Viburnum tinus, arbutus unedo, phillyrea latifolia, rubia peregrine, asparagus acutifolius, asplenium onopteris (148-149, Gamisans 1999) being just a few. These plants grow all over the island and make surveying on Corsica very difficult (16, Given 2004; 176, Jong et al 2001) because they are evergreen and occur in high densities.
Mediterranean ecological studies, whether directly monitoring the maquis or monitoring other events such as wildfires, always mention these scrublands or its plants (269, Lasaponara 2006). Remotely sensed studies of the maquis suggest that increased photosynthesis is related to increases in CO2 and precipitation rather than to factors such as climate change or seasonal variation (2276, Calvao and Palmerin 2011; 1904, Osbourne and Woodward 2001). These studies of the maquis also discuss the difficulty in classification of the maquis within remote sensing data. This arises because the maquis is comprised of multiple plants, each of which have unique spectral signatures, and because, unlike fields with defined borders, these scrublands taper in and out with varying density ( Jong et al 2001).
1.3 Normalized Difference Vegetation Index (NDVI)
NDVI is a representation of the amount of photosynthesis occurring within a particular pixel and depends heavily on the types of plants within the pixel (176, de Jong et al 2001; 269, Lasaponara 2006). This is because an NDVI can represent different physical things, such as one plant with many photosynthesizing leaves or many smaller plants. Deciding which is which is difficult from such an index and because of that NDVI often requires ground-proofing to make sure that the estimates made from those signatures accurately represents the realities of the area. NDVI uses the ratio between visible light and near infrared light to calculate an index for photosynthesizing plants. Chlorophyll absorbs visible red light for use in photosynthesis, so a healthy plant will reflect very little light in that range (Weier 1999). When seasons change, or other ecological factors cause plants to change their amount of photosynthesis the amount of reflected red visible light changes. However, the cell structure of leaves always reflects near-infrared light (Weier 1999). By subtracting the amount of red wavelength from the amount of infrared and dividing it by the addition of the two an index of photosynthesis is created. Places that have high photosynthesis, and are therefore absorbing most of the red wavelength, will have indices close to one since the reflected red wavelength will be so small. Places however that have decreased photosynthesis will have negative values because the red wavelength will be so much larger than the near-infrared wavelength.
Seasonal changes often cause the biggest differences between wavelengths but certain plants respond better to weather events, such as rain, and may experience a change in NDVI because of that. Tracking these changes may allow an experienced researcher to recognize plant species in an image if the researcher has knowledge of plants in the area. Other studies have suggested that these changes represent a change in land cover classification. While this may be true, the increase of NDVI in scrublands of the Mediterranean region related to an increase in precipitation and CO2 rather than land cover change (1904, Osborne and Woodward 2000). To use a combination of the NDVI and area represented by a single pixel may give some estimate of plant density.
1.3 MODIS These data are distributed by the Land Processes Distributed Active Archive Center (LP DAAC), located at the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center (lpdaac.usgs.gov). The Moderate Resolution Imaging Spectroradiometer (MODIS) is an instrument aboard two satellites the Terra (EOS AM) and the Aqua (EOS PM). The Terra satellite passes south to north across the equator in the morning and views the entire Earth‘s surface every one to two days. The Aqua satellite follows the same route but at night. The instrument collects 36 different bands, which monitor wavelengths from 0.4µm to 14.4 µm. The band numbers do not correlate with monitored wavelength, with some bands overlapping wavelengths (LP DAAC).
Bands 1-2 have up to a 250m spatial resolution, spectral range of 620-670 nm and 841-876 nm respectively; and are used for looking at land cover transformation, vegetation chlorophyll, cloud amount, and vegetation land cover transformation. Bands 3-7 have a spatial resolution of 500m, spectral range 459-479, 545-565, 1230-1250, 1628-1652, and 2105-2155 nm respectively; and are used to monitor soil and vegetation differences; green vegetation, leaf and canopy differences; snow and cloud differences; as well as cloud and land properties. In relation to archaeological research, MODIS‘s low spatial resolution makes it almost useless for most archaeological studies since anthropogenic features are rarely larger than a single MODIS pixel. However, this project follows in the footsteps of ecological research and the purpose of MODIS is to monitor ecological change (LP DAAC). Vegetation indices, such as MODIS product MOD13Q1, monitor the amount of photosynthesis occurring in a certain area. These changes can be responses to seasonal variation or other factors.
LP DAAC distributes processed MODIS Land Products that utilize MODIS raw data to create vegetation indices so that users do not need to process the raw data themselves. It also allows for consistency between processing so that studies using the same product, for this study the MOD13Q1, will be compare apples to apples. The MOD13Q1 is a vegetation index that uses “blue, red, and near-infrared reflectances, centered at 469-nanometers, 645-nanometers, and 858-nanometers, respectively” (LP DAAC) to determine the daily indices. This product comes with the NDVI raster, as well as an Enhanced Vegetation Index (EVI), a pixel reliability summary, and other data sets. These are provided every 16 days and use a pixel replacement method to eliminate cloud-covered images (LP DAAC).
2. Methodology and Workflow
The applicability of this method depends heavily on the availability of data for the study area. Previously, this method had hoped incorporate the Landsat-5 data into the workflow in order to get higher resolution data and perhaps incorporate a technique called spectral unmixing. However, neither USGS nor ESA, the two agencies that supply Landsat-5 data free of charge, were able to supply data for the study area and spectral unmixing when used on a MODIS pixel would uninformative since it is rare that a 250m pixel consists of a single endmember. Therefore, this project relied on MODIS data with a 250m pixel size and the ability to monitor changes in NDVI. This method requires the use of Marine Geospatial Ecology Tools (MGET), through Duke University, and accessible at http://code.env.duke.edu/projects/mget, and was executed using ArcGIS 10 (Roberts et al. 2010). MGET is required because MODIS data downloads as an .HDF file and ArcGIS does not read these files.
This project used 12 images, from December 2004- July 2010, obtained using the USGS Global Visualization viewer (GloVIS) or USGS Earth Explorer (EE) accessed through LP DAAC. A summer image, taken from July and a winter image, taken from November, December, or January, was obtained for any given year. Summer images were obtained for consistent dates, 193, but the increase in cloud cover during winter created variation in the selection of a winter image. These are listed in table 1 with the date given in a 365-day format, with “001” corresponding to 1 January. Once obtained, a header, describing what was in the .HDF file needed to be extracted using Marine Geospatial Ecology tool > Data Management > HDF Files > Extract Header in ArcGIS (Roberts et al. 2010). This header can be read through Notepad and data such as the name of the data set, the lower left corner, and cell-size need to be found or calculated from this data.
Since the same MODIS product, MOD13Q1 was used for every image in the workflow, the lower left corner and cell size for each data set was identical (LP DAAC). In the MOD13Q1 data these elements are labeled “westbound” for “X coordinate of lower-left corner” and “southbound” for “Y coordinate of lower-left corner”. The cell size was calculated according to the MGET website as (Northbound – Southbound)/ (the dim0: size) (Roberts et al. 2010). For this set of images the X coordinate was about 0, the Y coordinate was about 40, and the cell size calculated to about 2.08e-3. At this point additional information could be added in order to transform the data to a different coordinate system but that was not considered necessary. This data was then used to convert the SDS in .HDF to an ArcGIS raster.
There are at least two different types of rasters that can be extracted from the MODIS .HDF files, the “250m 16 days NDVI” and the “250m 16 days pixel reliability,” which are under the same title for each image. The NDVI raster will give the NDVI for each pixel and will have values that range from -3000 to 10000. However, valid values are only from -2000 to 100000 (LP DAAC). In order to get rid of the -3000 values one needs the pixel reliability raster. This raster, on a scale from -1 to 3, asses the reliable the data from each pixel. Anything with a value of -1 is not processed data and does not contain a calculated NDVI, values of 0 are good data, values of 1 are marginal data, values of 2 are snow or ice pixels, and values of 3 are cloudy pixels (LP DAAC). Particularly, pixels with a value of -1 or 3 will not be useful. Therefore, the reliability raster can be reclassified so that values of -1 and 3 will have NoData and values of 0-2 are 1.
Using the raster calculator, the reclassified pixel reliability raster and the NDVI raster may be multiplied together to eliminate water and unreliable pixels, those with a value of -3000, and then multiplied by the scale factor 0.0001 to get the values between -0.2 and 1, the valid range. Since we were interested in looking for pixels that went through significant changes in level of photosynthesis the scaled raster can be reclassified so that values between -0.2 and 0 are given a value of 1, 0- 0.66 are given a value of 2, and .66-1 are given a value of 3. Negative NDVI values are values where the reflected red wavelength is larger than the reflected near-infrared wavelength, which suggests that the pixel is not photosynthesizing very well. An NDVI of higher than .66 suggest that the majority of the red wavelength is absorbed and that the pixel is photosynthesizing very well. Values between 0-0.66 are less easy to define and a change from one of those values to either highly photosynthetic or scantily photosynthetic is not a drastic change.
After reclassifying all of the images to values of 1, 2, or 3, the raster calculator was used to subtract older images from younger images. It is important at this stage to make sure that the subtraction for each calculation is in sequential order but since we are interested in magnitude whether one subtracts older from younger or vice versa does not matter. This means that every raster, except for the December 2004 and the July 2010 rasters, will be used twice in a calculation, once as the younger image and once as the older image. This then produces a raster that represents the degree of change that occurred between one image and the next image. These values, in discrete integers, range from -2 to 2 with -2 and 2 representing pixels that went from highly photosynthesizing to scantily photosynthesizing or vice versa.
All of these images were then further reclassified so that values of extreme change, -2 or 2, were isolated from all other pixels. However, the coding for each image is different from the next. In order to monitor when in time extreme changes occurred, values of 2 or -2 were reclassified to a prime number unique to a single raster. All other pixels, ranging in value from -1 to 1, were given a value of 1 so that if change occurred at a later stage in the same spatial area it would not be lost. Table 2 lists the codes of prime numbers assigned to different rasters, the temporal changes, and their sequential number. After all rasters had been reclassified, all eleven rasters, representing changes from winter to summer and summer to winter from December 2004 to 2010, were multiplied together. This created a final raster with a set of unique values that could be broken down into prime number factors in order to examine when change occurred.
3. Analysis
3.1 Data Analysis
Of almost 5000 different possible values, the 54 that occurred are listed in table 3, along with what they mean in relation to chronology, and the number of pixels that exhibited that pattern. The majority of pixels, 99.4%, did not have any significant change occur in them. 842 pixels exhibited some sort of change and their unique patterns of change are represented by the remaining 53 values. 61.6% of these pixels changed only once but the single combination of change that occurred the most, 16.6% of the time, were significant changes between July 2008 to December 2008 and another significant change between December 2008 and July 2009. Four pixels exhibited a more advanced pattern, having had significant changes at least 5 times, and three of them exhibited change every round for the sequence 8-11. The time sequence where change occurred the most was between July 2007 and Januray 2007 or between December 2008 and July 2008. This was represented by the fact that the prime numbers 11 and 19 both had 18 values that contained them as factors.
Based on the data, it is clear that, according to our definition, very little seasonal change in photosynthesis has occurred in the last 5 years. Looking at figure 1 the majority of pixels that changed appear to be those along the coast or near where large patches of “NoData”, the white spaces, occurred. A few pixels do not follow this rule and are surrounded by pixels that had no change occur. These isolated patches are usually places where only one change occurred and may represent a burning event. The positioning of these pixels along the coast or near “NoData” regions can be explained by physical attributes. Examining the individual NDVI rasters there is significant variation in where the coastline and so variation in NDVI in these regions may be attributed to where the land “ended” in each image. There is also some evidence that the pixels near the coastline may be reflectances of large roads or cities based on the patterning of the pixels. Pixels near the “NoData” regions may be the effect of clouds that were not marked as such during the pixel reliability phase or where changes from snow pixels to photosynthesizing ones occurred.
3.2 Evaluation of the Method
One of the most frequent problems when executing this workflow was the automatic shutdown of ArcGIS 10 after some of the operations and the infamous 99999 error in which calculated rasters will not automatically display in ArcMap. This made the formation of a geoprocessing tool, which combined all of the above steps, useless since the model would usually fail after one process that involved a raster.
Another limitation of this method is that areas that had cloud covering in one image and none in another were excluded from study. Any raster calculation that involves a value of “NoData” just returns the “NoData” value and is not part of the calculation process. Since unreliable pixels, such as unprocessed or cloud covered, were given values of “NoData” based on the reliability raster, there were chunks of the study area that had “NoData” even though they could have, at least during some parts of the sequence, had the significant changes in NDVI. How to combine the rasters in a way that would include the reliable pixels available from areas that at one point have NoData, or classifying them in a way that would include only land “NoData” values, would be a significant improvement and may positively change our results.
Another limitation was the definition of a “significant” photosynthetic change. By classifying the data, we lost some information, such as the exact NDVI of a pixel. This choice meant that a pixel with an NDVI of .66 that switched to an NDVI of 0.0 would have the same “change” as a pixel with NDVI of .99 that switched to an NDVI of -0.2. It would also ignore when a pixel with an NDVI of .65 changed to an NDVI of -0.2, since the first NDVI would be classified as type 2 and then as type 1 and therefore an overall change of -1, even though its change in NDVI was of comparable magnitude as a switch from .66 to 0.0. It may be useful to, in the future, apply a different approach to classifying significant change, such as percent difference. This would allow for pixels that were otherwise classified as not having enough photosynthetic change, because they had a value between 0 and .66, to be recognized if they had a large enough percent change. Reclassification with prime numbers, perhaps differentiating between pixels with 60% change and not, could occur after this point.
It is also possible that we are not seeing much change in photosynthesis because of the coarse data set. 250m is a very large area and it is possible that there are patches within those blocks that have the sorts of cyclical photosynthesis we are looking for but they are small in comparison to the pixel. The use of data with a finer resolution may aid in differentiating those smaller areas.
4. Conclusion
This method of utilizing remote sensing does in fact locate places where there was significant change in NDVI. Areas that would have exhibited a pattern of change every round would have signaled areas that had specific high and low photosynthesis times, relating to a particular season, allowing archaeological survey during the low photosynthesis time, however, there were only three pixels that exhibited this sort of pattern. The majority of pixels fell into the expected category of “evergreen” probably because of the evergreen maquis that covers much of the study area. Combining this data with road and other logistical data would probably not result in a much more efficient survey plan since there are so few pixels that exhibited significant photosynthetic change. However, altering the method to use percent difference, exploring the NDVI of pixels that only changed once, incorporating “NoData” regions, and getting finer resolution data may be the next step in this procedure.
If these methods help to increase the amount of pixels with significant change in NDVI or help us to isolate what specific changes mean this method of planning targeted survey regions would be extremely useful. One could also combine this data with the most recent MODIS NDVI product to isolate and survey pixels that have low NDVI at the time during the study time. After making these changes, transforming the MODIS data into a more widely used map projection, such as UTM, and combining it with logistical data, such as roads and cities, would allow for a comprehensive “survey plan” based on NDVI to be put in place.
References
Calvão, T. and J. M Palmeirim. 2011. A comparative evaluation of spectral vegetation indices
for the estimation of biophysical characteristics of Mediterranean semi-deciduous shrub
communities in International Journal of Remote Sensing, 32(8): 2275 -2296.
De Jong, S.M, T. Hornstra, and H.G. Maas. 2001. An integrated spatial and spectral approach to
the classification of Mediterranean land cover types: The SSC method in International
Journal of Applied Earth Observation and Geoinformation, 2: 176-183
Fentress, E. 2000. What are we counting for? in Extracting Meaning from Ploughsoil
Assemblages (eds R. Francovich, H. Patterson, and G. Barker). Oxford: Oxbow, pp.44
-52.
Gamisans, Jacques. 1999. La Végétation de la Corse. Aix-en-Provence: Édisud, 143-153.
Given, M. 2004. Mapping and manuring: can we compare sherd density figures? in Side-by-Side
Survey: Comparative Regional Studies in the Mediterranean World (eds. S.E. Alcock and
J. F Cherry). Oxford: Oxbow Books, pp. 13-21.
Kuemmerle, Tobias , Röder, Achim and Hill, Joachim. 2006. Separating grassland and shrub
vegetation by multidate pixel-adaptive spectral mixture analysis in International Journal of Remote Sensing, 27(15): 3251-3271.
Kvamme, Keith. 2006. Magnetometry: Nature‘s Gift to Archaeology in Remote Sensing in
Archaeology: An Explicitly North American Perspective, Tuscaloosa : University of Alabama Press, pp. 205-235.
Lasaponara, R. 2006. Estimating spectral separability of satellite derived parameters for burned
areas mapping in the Calabria region by using SPOT-Vegetation data in Ecological
Modeling, 196: 265-270.
Llobera, Marcos, K. M. Wilkinson, M. C. Weiss, R. J. Flaming, N. A.F. Marini, and S. Mazet.
2010. Into the maquis: Methodological and interpretational challenges in surveying la
Balagne, northwest Corsica in Journal of Mediterranean Archaeology, 23(2): 169-196.
Marshall, A. 1998. Visualising burnt areas: patterns of magnetic susceptibility at Guiting Power
1 Round Barrow (Glos., UK) in Archaeological Prospection, 5: 159-177.
Melendez-Pastor, Ignacio Melendez-Pastor, Ignacio, Jose Navarro-Pedreño, Magaly Koch, and
Ignacio Gómez. 2010. Multi-resolution and temporal characterization of land-use
classes in a Mediterranean wetland with land-cover fractions in International Journal of
Remote Sensing, 31(20): 5365-5389.
NASA Land Processes Distributed Active Archive Center (LP DAAC). MODIS Vegetation
Indices 16-Day L3 Global 250m: MOD13Q1 . USGS/Earth Resources Observation and
Science (EROS) Center, Sioux Falls, South Dakota. 2004-2011.
Osborne, C.P. and F.I. Woodward. 2001. Biological mechanisms underlying recent increases in
the NDVI of Mediterranean shrublands in International Journal of Remote Sensing, 22(10): 1895-1907.
Roberts, J.J., B.D. Best, D.C. Dunn, E.A. Treml, and P.N. Halpin. 2010. Marine Geospatial
Ecology Tools: An integrated framework for ecological geoprocessing with ArcGIS, Python, R, MATLAB, and C++ in Environmental Modelling & Software, 25: 1197-1207.
Saturno, William, Thomas Sever, Daniel Irwin, Burgess Howell and Thomas Garrison. 2007.
Putting us on the Map: Remote Sensing Investigation of the Ancient Maya Landscape in Remote Sensing in Archaeology (eds James Wiseman and Farouk El-Baz). New York: Springer, pp. 137-160.
Slater, L.A., N.D Hamilton, A. Sandberg and M. Jankowski. 2000. Magnetic prospecting at a
prehistoric settlement in Main in Archaeological Prospection, 7: 31-41.
Weier, John and David Herring. 1999. Measuring Vegetation (NDVI & EVI) on NASA Earth
Observatory, web. 9 Dec 2011.
