GIS 5990 Final Project: Predictive Model of Neolithic Sites in Southeast Turkey

Figure 1: Overview of study area.

The aim of this project was to create a predictive model of Neolithic archaeological sites in southeastern Turkey. This area was chosen due to the location of Gobekli Tepe in the region. This site has become well known as one of the earliest uncovered monumental structures with some researchers (such as excavator of the site Klaus Schmidt) hypothesizing that mobile hunter gatherers constructed the site. The impetus for its construction may have been a change in ideology that also created fertile ground for the transition to sedentary, agricultural lifeways.

Figure 2: Site in the study area.

A successful predictive model would assist in explaining site location and help researchers in discovering new sites. Site data were acquired from The Archaeological Settlements of Turkey GIS (http://tayproject.org/giseng.html). Neolithic site data were manually entered into a new shapefile in ArcMap. This data is the weakest link in the analysis due to the potential for spatial inaccuracy.

The digital elevation model (DEM) seen in figures 1 and 2 were created by using the Mosaic tool on fifteen ASTER images acquired from the USGS Earth Explorer website (http://earthexplorer.usgs.gov/). With the image properly projected as a DEM, slope and aspect data could be generated in addition to elevation data.

Figure 3: Elevation data reclassified.

Based on the location location of Neolithic sites at lower elevations, the Reclassify tool was used to classify the elevation data into four equal categories. Category 4, representing the category with the highest probability of sites, was from 0-800 meters. In descending probability, category 3 was 800-1600 meters, category 2 was 1600-2400 meters, and category 1 was 2400-3200 meters.

Figure 4: Original slope data.

Figure 5: Original aspect data.

Figure 6: Slope data reclassified.

Figure 7: Aspect data reclassified.

Slope and aspect data were similarly classified into four categories of high to low site probability. The highest probability category for slope was defined as 0-5%. In descending order of probability, the remaining three categories were 5-15%, 15-25%, and 25-75%.

An aspect of 0 (or flat) was rated as the highest probability, with southern-facing aspects being the next probability category. East and west facing aspects were the next category, with northern facing aspects being of the lowest probability.

Figure 8: Water buffer of 500 meters.

Data on rivers and lakes, acquired from Eurostat (http://epp.eurostat.ec.europa.eu), were buffered to 500 meters. A raster was then created by using the Feature to Raster tool with probability defined as binary categories. Areas within the 500 meter buffer were considered high probability with all else being low probability.

Figure 9: Weighted Overlay.

After the data on the four variables were suitably processed, the Weighted Overlay tool could be used. Elevation and water source distance were each weighted 35% while aspect and slope were each weighted 15%. The result can be seen above.

Figure 10: Ordinary Least Squares analysis.

In order to test the weighted overlay predictive model, an ordinary least squares analysis would be conducted. First, the Create Random Points tool was used to created over 100 points to function as locations without sites. This was merged with the original site shapefile. The Extract Multi Values to Points tool was used to incorporate data from the variable rasters into the merged point shapefile. The Ordinary Least Squares tool can then be run.

Unfortunately, the statistics indicate the model is not successful in explaining site locations. Only elevation returned as a significant variable. There are several explanations for why the other variables were not significant. First, as noted above, site location accuracy is not assured based on the data source used. A more robust model would incorporate verified site locations. Additionally, the reclassifications of slope, aspect, and water source distance may have been incorrect. Other classifications may function better. It is also possible that slope and aspect do not in fact help explain site location.

We might also attempt to incorporate paleoenvironment data, including on water sources. Using modern hydrography as a proxy for ancient water sources may not be justified.

Another potential variable to include in future models is viewshed and site visibility. Figure 11 shows an experiment with the Observation Points tool depicting raster cells visible from five well-known sights in the study area. Gobekli Tepe in particular has a wide view to the south. There is potential here for future model attempts.

Figure 11: Observation points for five site in study area.

GIS 4035 Module 10: Supervised Classification

Supervised classification and spectral distance file of Germantown, Maryland.

The task for module 10 was to utilize ERDAS Imagine to create, analyze, and edit supervised classifications of multispectral images. We began by reviewing the fundamentals of the Signature Editor tool (to create and edit the training samples used for a supervised classification) and creating an AOI (area of interest) layer. We then used the Inquire tool to locate specific points on the image on which we created polygons to capture to areas to be used as the spectral signature of a class. We also reviewed the use of the Seed tool; this tool automatically creates a polygon capturing pixels of a set spectral distance from the initial point based on input parameters. I preferred using the Seed tool to manually inserting polygons as it gives more control over keeping each class distinct. We also spent some time analyzing the histograms and mean plots of spectral signatures. These tools give use a more precise way to ensure our classifications capture unique features and are not spectrally confused.

The assignment required us to use all the tools reviewed in the module to create our own supervised classification of an image of Germantown, Maryland. We were given coordinates of three urban features, two fallow field features, four agriculture features, and one each for grasses, deciduous forests, and mixed forests. We also were to create signatures for water and roads. With this completed, we then needed to recode the number of classes down to eight, consolidating the categories with multiple classes into one each. We also needed to include the spectral distance file in our final map; this image shows us at a glance the pixels in the image that are furthest away from any of our classes. The final layout was completed in ArcMap.

GIS 5990 Module 9: Biscayne Shipwrecks, Analyze Week

Five shipwrecks in Biscayne NP with 300 meter buffers showing benthic zones.

Reclassified layers showing benthic zone and bathymetric data; these layers were used in the weighted overly below.

The output of the Weighted Overlay tool using the reclassified benthic zone and bathymetry layers.

This was the second week of our three-week project focusing on GIS use in analyzing shipwrecks in Biscayne National Park. This module was similar to our previous work in creating predictive models, only this time we processed data relating to the sea floor. As seen above, the variables included in our weighted overlay were benthic zones and bathymetry. Interestingly, the bathymetric data seemed to contradict the actual locations of known shipwrecks; most shipwrecks were not located in the shallowest waters. To account for this, the weighting of the bathymetric data was lowered to 30%; the benthic zone data comprised the majority 70% weighting.

For the most part, this was an enjoyable assignment to complete. I did run into a frustrating difficulty at the end as I attempted to run the Weighted Overlay tool. Many, many attempts either ran into errors or resulted in a raster output that was not classified correctly. Eventually, I checked the properties of my input data, the reclassified benthic zone and bathymetry rasters. I discovered the latter raster was lacking a defined coordinate system; fixing this allowed the tool to run quickly and correctly. The experience reinforced the lesson that I must always check my data to ensure everything is in order.

GIS 4035 Module 9: Unsupervised Classification

Five category unsupervised classification of the UWF campus based on true color imagery.

The ninth module of the course guided us through performing unsupervised classifications of imagery in ArcMap and ERDAS Imagine. An unsupervised classification takes basic guidelines from the user (such as the number of desired categories) and creates categories based on the appearance of each pixel. A perfect classification would, for example, classify all water in a category, all trees in another category, and so on. There has likely never been a perfect classification, however, so the resulting classified image must be edited to better capture the desired categories.

In order to create the above classified image of the UWF campus, the original true color image was run through Imagine's Unsupervised Classification tool to create fifty categories. The resulting image (not shown) looked very similar to the original image. We then reclassified each of the fifty categories into the five classes of trees, grass, buildings/roads, shadows, and mixed (grass/urban). The main source of error was the overlap of bright grass and ground areas to some urban areas; this created the need for the "mixed" class. 

GIS 5990 Module 8: Biscayne Shipwrecks, Prepare Week

Five Biscayne National Park shipwrecks in relation to nautical charts (modern and historic) and bathymetry.

Module 8 begins our three week investigation of shipwrecks within Biscayne National Park. This week of preparation involved locating, downloading, and properly preparing bathymetric data as well as historic nautical charts. The National Oceanic and Atmospheric Administration was the source for both. The Historic Map and Chart Collection (http://historicalcharts.noaa.gov/) provides a wealth of historic charts, with the only caveat being they must be georeferenced once acquired. Although initially frustrating, this proved to be simpler than I had assumed. Having the data represented in the historic chart (in this case from 1856) is important if one is interested in knowing the quality of data available at the time of a shipwreck. 

The bathymetric data proved more troublesome to locate. After trying several downloads (including multibeam data that I may not have processed properly), I finally chose the image above featuring good detail but covering only the western portion of the Park. This image as well as the historic chart were clipped to the boundary of Biscayne National Park. 

This was an enjoyable module that proved slightly frustrating when it came to finding the best data source. Nonetheless, I am glad to get more experience in locating and processing data for a project.

GIS 4035 Module 8: Thermal and Multispectral Analysis

Imagery combining the thermal infrared band (6) with the shortwave infrared bands (5 and 7) to highlight fires.

The eighth module of the course focused on using thermal infrared imagery to extract data unavailable from other EMR wavelengths. The general concepts and tools were similar to previous exercises, yet thermal infrared imagery presents unique challenges. First, thermal infrared EMR is emitted, not reflected; the amount of thermal infrared EMR emitted by a feature is a combination of the amount of energy absorbed, the feature's composition and surface characteristics, and the sensitivity and exposure length of the image (among other factors). In simple terms, thermal infrared EMR reflects temperature, but we cannot assume a direct correlation without calibration.

Our deliverable for the week was more open-ended than usual; we were to choose an area or feature in one of the two composite images created for the module and create an image that highlights the chosen area or feature. The thermal infrared band needed to at least be used to identify the feature even if it wasn't used in the final layout. While we had already identified select fires in the imagery, I was struck by how defined the fires' core extents appeared when the thermal infrared band was combined with the shortwave infrared bands. After applying a Gaussian stretch in Imagine and a minimum-maximum stretch in ArcMap, the fires popped out from the background significantly. This did allow me to identify a third, small fire southeast of the largest fire that I had not noticed previously.

I enjoyed the experimental nature of this module and learning of the unique data one may extract from thermal infrared imagery. I am still not confident in my ability to manually manipulate histogram breakpoints to achieve fruitful results, but this module did help me significantly understand the areas I need to work on.

GIS 5990 Module 7: Scythian Mounds Report Week, Regression Analysis

Results of ordinary least squares regression analysis on Scythian mound data.

The third and final week of our work with Scythian mounds near Tuekta, Russia allowed us to dig deeper into statistical analyses. I am glad to get more practice with statistics; perhaps eventually I will come to understand statistics and how to apply statistical tests in practical situations. Until then, the detailed guidance in our assignment instructions on interpreting statistical test results was most welcome.

We began the assignment by creating 100 random points as examples of locations without sites; these were then merged with our point shapefile digitized previously representing mound locations. We then needed to include the category of elevation, slope, and aspect in which each point was located in the merged shapefile's attribute table. I used the Extract Multi Values to Points tool, as described in the instructions. This was much simpler than manually entering the data, although it did entail deleting superfluous fields in the attribute table.

After the data were prepared properly, we ran the statistical tools of Ordinary Least Squares Regression, Spatial Autocorrelation, and Hot Spot Analysis. The results of these analyses indicate all three environmental variables affect mound placement (with elevation being the most significant) and give a 99% confidence value that mound placement is not random. The detailed results of the Spatial Autocorrelation analysis is included in the above map.

This module taught me much on the proper way to set up regression analysis. I will likely use this module for reference in my internship and future work until I become more comfortable with statistical work.

GIS 4035 Module 7: Multispectral Analysis

Feature 1 multispectral analysis: Deep water.

Feature 2 multispectral analysis: Snow.

Feature 3 multispectral analysis: Shallow water.

Our task this week was to practice using different methods of multispectral data analysis to find three features with particular EMR signatures. The methods practiced in the assignment included histogram analysis and manipulation, displaying different combinations of bands as composites, creating indices based on different bands (e.g., Normalized Differential Vegetation Index), and using the Inquire Cursor to get detailed data on individual pixels. In order the find the three features, I primarily used the Inquire Cursor method based on histogram analysis.

The assignment helped me greatly in understanding histograms, although I am still not confident in manipulating them for a better image. By looking at the histograms of each band, one gets a good sense of how the data will display. I also appreciated the practice in combining bands into custom composites to highlight particular data. This was something I did not quite grasp in ArcMap, but the methods in Imagine made the concepts clearer.

GIS 5945: Dream Job

It is difficult for me to pin down exactly what my dream job would be.  Any archaeology position in the western portion of the United States (and I even include Alaska) could be a dream job. However, I am not limiting myself to archaeology positions.

A job I recently applied for on USAJobs.gov was with the Forest Service in the Gila National Forest, Silver City New Mexico. The official job title is "Interdisciplinary Soc/Bio/Phy Information Specialist." The job description is to support "managers and resource specialists in the use of Geospatial Technology and Geographic Information Systems (GIS) for land, natural resources and ecosystem management." It is the interdisciplinary aspect of the job that appealled to me; working with a variety of experts from different disciplines to help solve problems and overcome challenges using GIS. It also helps that I am familiar with the region; the archaeology field school I attended for two summers was half an hour from Silver City. It is a region I came to enjoy and would not mind at all living in.

Because of its interdisciplinary nature, the job application requirement list consisted of numerous "or" statements. Essentially, the job required a Bachelors degree in a hard science, resource managment, or a social science (or equivalent experience). It further required one year of graduate education, specialized experience, or "superior academic achievement." I would qualify under "superior academic achievement," and my graduate coursework plus my internship would give me some additional points. 

The archaeology-focused courses in the UWF program would likely be the most beneficial to performing this job. While the projects may be more natural resource oriented rather than cultural resource, many of the same tools and techniques would apply. Although remote sensing expertise was not mentioned specifically, such a skill might also prove useful in the position.

GIS 5990 Module 6: Scythian Landscapes, Analyze Week

Geographic variables and Scythian mound locations.

This module continued our exploration of Scythian burial mounds near Tuekta, Russia. The majority of the assignment involved creating layers representing geographic variables that may have affected where mounds and other sites were placed. As seen in the map above, these variables include the topography, slope, aspect, and elevation of the region. The raw rasters representing slope, aspect, and elevation were reclassified to show the regions most and least likely to feature sites. The most favorable areas will be relatively flat, low elevation, and facing south (if they are on a slope at all). We also digitized points representing the centers of mounds in the provided aerial imagery. I was surprised to find so many, as previously I had focused mostly on the line of largest mounds in the eastern portion of the image. I ended up digitizing 66 mounds, although there are likely more. The pattern immediately seen is that the mounds form a series of lines running roughly north-south, with most of the lines trending towards a NNE-SSW orientation.

GIS 4035 Module 6: Spatial Enhancement

Landsat 7 image run through Fourier transformations, sharpening, and statistical filter in ERDAS Imagine.

Past modules in this and other courses have mentioned that GIS operators often must correct acquired imagery for various errors prior to using it for their analyses. This module introduced us to some of the ways such corrections can be done. We explored the tools available in both ERDAS Imagine and ArcMap, including Fourier transformations, high pass filters and low pass filters. Low pass filters generate output that appear smoother and less detailed than the original image; noise in the image is also removed, according to the size of the kernel chosen. High pass filters create high contrast and noisy output that highlights the edges of features.

For the assignment, we were to experiment with Imagine and ArcMap to create an enhanced version of a Landsat 7 image that minimizes as much as possible the scan line corrector failure striping without removing too much detail. The assignment walked us through using the Fourier Transform Editor tool as the initial step in the process. While I experimented with different wedge placements in the editor, I was unable to create an image that was significantly better than my first attempt. After experimenting with various filters and settings in Imagine and ArcMap, the final image used above was the result of sharpening and a statistical filter in Imagine. Other filter and setting combinations created images similar to this, and a few were significantly worse. I did not experiment much with the histogram, however; I hope to learn more about histogram manipulation in future modules.

GIS 5990 Module 5: Scythian Landscapes, Prepare Week

Layout from the first week of the Scythian Landscapes project.

After our two week module on predictive modeling, we begin another three chapter project investigating a particular topic. This time, we are considering the theoretical perspectives contained within landscape archaeology and, with GIS as the tool, looking at Sythian burial mounds in southern Russia (near the borders of Mongolia and Kazakhstan).

Our tasks for the first week were relatively light. To prepare for next week's tasks, we prepared a mosaic of four DEMs that were then clipped to a polygon representing the study area. We were provided with an aerial photo of the mounds near Tuekta, Russia that needed to be georeferenced to the clipped DEM. Our layout was required to include the clipped DEM along with an inset of the georeferenced aerial photo. I added a small-scale reference inset to further orient the viewer; the arrows I added from Adobe Illustrator to make the relationship between the insets clear.

GIS 5990 Module 4: Predictive Modeling

Weighted overlay depicting a predictive model for site locations.

Module 4, stretching over two weeks, required us to progress through the multi-step process of creating a predictive model for archaeological site location probability. A predictive model combines known site locations with variables (generally environmental) to create a classified raster categorizing the probability of finding sites in different regions. Common variables to consider are elevation, slope, aspect, and distance to water sources; these were the variables used for our project. In order to create data layers that can be used in the Weighted Overlay tool to produce the predictive model, the required data must be located, acquired, clipped, merged, reclassifed, and/or converted. As an example, we must decide how many classifications we will have in our weighted overlay layer; we used three for our assignment. Thus, once we extracted an aspect layer from the DEM using the Aspect tool, we then used the Reclassify tool to create a new aspect raster with only three classifications. We had to similarly edit and create the slope, elevation, and water source distance layers for use in the Weighted Overlay tool.

This assignment, although lengthy, was an enjoyable practice of something I have been doing in my internship with the BLM. Something both this project and my internship have driven home is the need for a good project plan and good data for any predictive model. The sheer number of variables and different potential classifications and weightings for those variables, as well as varying responses to those variables by different human populations in different times and environments, requires careful consideration if one has any hope of producing a useful predictive model.

GIS 4035 Module 5a: Introduction to ERDAS Imagine and Digital Data

Subset of a classified image selected from ERDAS Imagine and exported to ArcMap for layout finalization.

Our module this week was divided into two topics. First, we dove deeper into the details of the electromagnetic spectrum, including the relationship between wavelength, frequency and energy along with how to calculate each. While we need not be physicists to use remotely sensed data, we do need to have an understanding of the structure of the electromagnetic spectrum and the type of data we can extract from different wavelengths.

Our second topic was an introduction to the ERDAS Imagine program. I am glad to finally be learning this program, in part due to the positive things I've heard of it from others but also as a change of "scenery" from ArcMap. Our deliverable (above) was simply a selected section from a classified image exported and finalized in ArcMap. Most of the assignment was designed to get us familiar with the program for use throughout the semester and, hopefully, our careers.

An issue that often arose as I worked with ERDAS Imagine was the difference between Imagine 2011 (the basis for the assignment) and Imagine 2014. The majority of changes I encountered were simple to recognize and overcome. The only major problem I had was in exercise 3; the subset image I created using the inquire box never updated the area attribute. I made several attempts without success. Most likely I was in error but I could never locate the source. In the end, I calculated the updated area totals in ArcMap using the Count field of the attribute table.

I enjoyed working with Imagine as the program seems to be much more efficient than ArcMap in simple manipulations of raster data (such as panning and zooming). I am looking forward to learning more about what the program can do.

GIS 4035 Module 4: Accuracy and Ground Truthing

Land use land cover classification of Pascagoula, MS with accuracy of thirty sample points symbolized.

Our task for week 4 of Aerial Photo Interpretation and Remote Sensing was to check the accuracy of last week's efforts to classify an aerial photograph. The best methods for ground truthing such a classification involve direct testing sample sites in the field, but the online nature of the course prevents such tests. Instead, we utilized Google Maps (and especially Street View) to test the classification accuracy of thirty sample sites. I roughly followed a stratified random sampling pattern based on classification category. However, small or homogeneous categories were given fewer samples (e.g., the large body of water in the west, the cemetery) while large, heterogeneous categories were allotted more samples (e.g., residential, commercial and services). The accuracy of my classification turned out to be about 73%. The main error sources were misclassifications of bodies of water and of forest cover. While my misclassification of Krebs Lake as a bay could, perhaps, be forgiven, my misclassifications of forest cover were due to not correctly distinguishing deciduous and evergreen trees. As a result of this assignment, however, I believe I would have a higher accuracy percentage on a similar aerial photo.

GIS 5990 Module 3: Mayan Pyramids and Angkor Wat, Analysis and Report

Supervised classification of Angkor Wat, Cambodia based on SWIR composite.

Google Earth image of potential pyramid locations in El Mirador.  Potential sites, difficult as they are to see, are in red.

Google Earth image of various layers of El Mirador overlaying the standard Earth view.  The NDVI layer is visible.

Google Earth image of El Mirador supervised classification overlay with legend.

The final module of our initial project involved two main tasks.  First, we were to export the layers created in the first two modules into a format that could be displayed and shared in Google Earth.  This was a relatively simple task using the 'Layer to KML (Conversion)' tool for individual layers and the 'Map to KML (Conversion)' tool for entire data frames.  The results can be seen above.  I was disappointed in the appearance of my potential pyramid layer and some of the colors of my classification layer are off, but in general the process was a success.

A second task required of graduate students was to repeat the supervised classification process with Angkor Wat, a well-known archaeological complex in Cambodia.  We were to locate and download appropriate Landsat imagery, then classify the image based on a false color, NDVI, or SWIR image.  As La Danta pyramid was the focus of our El Mirador classification, the monumental core of Angkor Wat would be our focus here.  The difficulty resided in the different land cover types that are included within the core.  Every training sample I tried over-represented the likely locations of additional remains.  However, my best classification restricted the core sample to the vegetation overlying the structures.

This was an enjoyable module that would have been more so if my computer had not expired in the middle of the assignment.  Even so, I believe I was mostly successful in completing the two tasks.  As always, however, more time could be spent on tweaking the training sample for a more refined classification.

GIS 4035 Module 3: Land Use Land Cover Classification

Aerial photo of Pascagoula, MS with land use/land cover classification overlay.

This week for our Photo Interpretation and Remote Sensing lab we were given the challenge of classifying an aerial photograph based on land use and land cover.  While we were not required to get too detailed with our classification, we were required to classify everything.  We were to create a new shapefile and then create polygons over each classification.  I extensively used the "trace" and "clip" tools from the editor toolbar as I classified the image to ensure everything was properly classified with no gaps remaining between polygons.

This was a rather challenging lab assignment, due both to its intrinsic difficulty as well as to the catastrophic failure of my computer.  It was poor timing, but I believe I was successful with the above map in the end.

GIS 5990 Module 2: Identifying Mayan Pyramids, Data Analysis

The second week of our GIS Special Topics in Archaeology course continues our experimentation with remotely sensed data from Landsat 7 in our search for Mayan pyramids.  In this assignment, however, we focused on the infrared bands in looking at the electromagnetic radiation signature of the region's vegetation in the hope of locating archaeological remains hidden underneath.  We first utilized our false color composite from last week to create a normalized difference vegetation index (NDVI), which represents the difference between the red and infrared bands.  This is intended to show the plant life as communicated by EMR.  The next map we created was a composite consisting of Landsat bands 4,5, and 1.  Incorporating two bands of infrared data, this composite represents more detailed information on the health of plant life.  The final map is the result of a supervised classification of the 451 composite, created from a training sample focused on locating more remains based on the EMR signature of the La Danta pyramid.

I made several goes at the training sample before I reached one I was reasonably satisfied with.  My final classification over-represents bare ground, but it is less egregious than my earlier ones.  Not many areas returned positive hits for pyramids which is better than too many hits.  I do believe the areas returned as possible archaeological remains would be worth a look, although further data and analysis would always be welcome.

GIS 4035 Module 2: Visual Interpretation

The first two modules of our Photo Interpretation and Remote Sensing course focused on the background and fundamentals we need to build on throughout the semester.  In particular, the two maps above show our experimentation in module two with the different criteria by which we interpret remotely sensed data (specifically aerial photography).  In the first map, we were required to identify five categories of tone (from very light to very dark) and five categories of texture (from very fine to very coarse).  Our next task, reflected in the second map, was to use four additional identification criteria (shape and size, shadow, pattern, and association) to interpret the photograph.  The most difficult aspect of this task was attempting to consider each criterion in isolation; in normal interpretive practice we obviously use all available criteria.  However, considering them singly was helpful; I had not considered how useful shadows might be to correctly interpreting aerial photography.

These introductory modules provided a good foundation in the science of remote sensing and the basics of interpreting the resulting data.  I look forward to learning much more throughout the semester.

GIS 5990 Module 1: Identifying Mayan Pyramids, Data Preparation

The first module of our Special Topics in Archaeology course begins a three module project using remotely sensed data to locate Mayan pyramids in dense jungle terrain.  As seen above, we produced three rasters depicting a single area in El Mirador, Guatemala, with the intent of noting any noticeable characteristics in the raster data indicating the presence of the La Danta pyramid.  The Landsat 7 data were downloaded from the USGS Earth Explorer website (http://earthexplorer.usgs.gov).  The first raster we looked at was composed solely of Landsat 7's band 8.  This panchromatic (sensitive to visible light) band is of a higher resolution than the other bands (15 meter versus 30 meter), yet the La Danta pyramid is far from obvious in the raster.  We next combined Landsat 7 bands 1, 2 and 3 into a composite raster to replicate natural colors.  Again, the pyramid does not show up as terribly obvious.  The finally raster produced this module was another 3 band composite, only this one utilized bands 2, 3 and 4 to create a false color raster showing near-infrared data (band 4) as red in the final image.  While perhaps a slight improvement over the previous two rasters, I was still unable to clearly identify the pyramid in the false color raster.

In addition to the bands utilized in this module, Landsat 7 data include three additional bands of short-wave and thermal infrared remotely sensed data.  Further raster composites are possible utilizing these bands for more analytical tasks in the infrared wavelengths.  Perhaps these bands will draw out the La Danta pyramid more clearly in the next two modules.

GIS 5103 Module 11: Sharing Tools

The eleventh and final module of the course taught us the correct procedure to create easier to share tools.  The ability to share script tools with others is one of their key benefits, yet it is important to provide necessary documentation as well as the correct file structure to ensure it will run properly.

For our assignment, we were provided with a script, script tool and toolbox.  While the script tool parameters were preset, we were required to edit the script itself to accept user inputs for two of the parameters.  Because we were also instructed to rename the provided files, the script tool properties needed to be updated in order to locate the script.  We also needed to edit the script tool's item description for each parameter in ArcCatalog.  This significantly improves the user interface for those not familiar with the tool, as the description we entered would be displayed in the help window of the script tool's dialog box.

Finally, we embedded the script into the tool and set up a password.  This allows the tool to be shared more easily by dispensing with the need for the stand-alone script.  It also increases security on the tool by preventing unauthorized users from altering the script.

GIS 5103 Participation Post #2: Viewshed Analysis in Archaeological Modeling

In their 2014 article "A house with a view? Multi-modal inference, visibility fields, and point process analysis of a Bronze Age settlement on Leskernick Hill (Cornwall, UK)" (2014, Journal of Archaeological Science 43, 267-277), Stuart Eve and Enrico Crema describe their efforts to model Bronze Age sites using different sets of variables and critique the way in which some researchers approach statistical modelling.  Modelling, as the authors use the term, involves performing spatial analyses on variables that potentially affected site placement; a successful model would reveal correlations between variable values and site placement.  For example, a model might propose that sites are more likely to be located at certain elevations or in areas of a certain level of precipitation.  Environmental variables such as elevation, slope, aspect, land cover, rainfall, and distance to permanent water sources are all common variables to consider when constructing a model for analysis.  Eve and Crema emphasize that no model is "true"; there are only more or less successful models.  Even a successful model, showing strong statistical correlation between variables and site placement, is not guaranteed to be closer to true than other models.  A wide range of variables and different levels of analyses must be considered.

Eve and Crema investigate three proposed models for the location of sites on Leskernick Hill.  The two models unique to their study involve line of sight analyses.  The first model proposed the placement of sites to maintain line of sight with ritual sites in the area while the second proposed a line of sight preference for tin deposits.  The third model looked at the more standard topographic and environmental variables.  The line of sight analyses were created by calculating the viewshed of every raster cell coinciding with Leskernick Hill using GRASS GIS and Python for batch processing.

The results indicate the need to consider difference scales of analyses.  While the first model (line of sight with ritual sites) was the best match overall, sites in the western portion of the study area fit the second model (line of sight to tin resource) much better.  Removing these western sites also showed the southern sites fit the first model even better.  The topographic inputs of the third model did not fit site location very well.

Link to article (hopefully one that works this time):

10.1016/j.jas.2013.12.019

GIS 5265 Module 10: Final Project

For our final project of the course, we were to either choose an independent project or work with the Oaxaca Valley data from module 6 for further catchment analyses.  As I could not find a suitable topic (or at least one for which data were readily available), I went with the Oaxaca Valley data.  Looking through the data provided for the northern portion of the valley (grid rows N11 through N14), there appeared to be oscillations in population density from the south to the north and back again through the time periods.  My project was to calculate maize consumption estimates based on population and compare this with estimated maize yields based on 1 km catchments.  

To create the Thiessen polygons as described in the module video for catchment analysis, I set up a model in ModelBuilder to create the final 1 km catchment buffer using the Buffer, Feature to Point, Thiessen Polygons, and Intersect tools.  I also used the model to create additional buffers of 2 km and 5 km, although I did not end up making use of them for this project.

The population data for my collection units for each period was acquired by joining the Excel sheets provided in our module 6 data.  Using the high and low consumption estimates provided in our readings (.16 metric tons and .29 metric tons respectively), I added four fields to the attribute tables for each period and used the Field Calculator to compute consumption estimates for each collection unit.

Arriving at maize yield estimates for the 1 km catchments was also a multi-step process.  The Identity tool was used to create shapefiles for each land type corresponding to the catchment areas.  A new field as added to calculate the area in hectares (since the estimates used in our readings are based on hectares).  Based upon data provided by Anne Kirkby (1973, The use of land and water resources in the past and present Valley of Oaxaca, Mexico, University of Michigan Museum of Anthropology Memoir 5), the yields of maize through time were multiplied by the area of each land type (again using the Field Calculator) to arrive at potential yields.

Microsoft Excel's Data Analysis correlation tool was used to determine if there was a correlation between population consumption estimates and maize yield estimates.  Unfortunately, no such correlation was seen.  The resulting data does show an oscillation of population from south to north and back, but significant variables remain unknown.  We do not know the fallowing strategy that was applied, nor do we know the potential degradation in catchment yields as the population grew and cultivation was intensified.  We also do not know the details of rainfall variation.  The majority of the study area (excepting one grid square) is considered to feature low average rainfall (below 700 mm), and the potential yields of the land types in the region (Type II and Type III) are highly susceptible to precipitation variation.  While this is suggestive as a potential cause for population movement through time, the lack of more detailed rainfall data prevents me from offering more than a mere suggestion.

GIS 5103 Module 10: Creating Custom Tools

Script tool messages after replacing print with arcpy.AddMessage().

Dialog box of script tool parameters.

While Python scripting for ArcMap allows one much more freedom to create functions not otherwise available, using and editing Python script requires a degree of knowledge and experience.  Fortunately, it is possible to import a stand-alone Python script into a toolbox as a script tool that is more easily shared with and used by others, regardless of their knowledge of Python.  Our exercise and assignment this week taught us to create and edit such script tools.

Although creating a script tool is a relatively simple matter, it does require a few steps to ensure the tool works properly without further editing the script itself.  After importing the script as a script tool into a toolbox, one must set up the properties and parameters of the tool.  The parameters include all the variables for the tool (inputs, outputs, etc.) as well as their data types.  This sets up the dialog box within the ArcMap GUI, but the Python code itself must be edited to accord with these parameters to allow users to successful enter variables into the tool.  This is accomplished by using the GetParameter or GetParameterAsText functions.  These are numbered in parentheses in the order they appear in the dialog box, starting with zero.  This will allow the script tool to function with custom inputs from the user.  The only remaining task is to use the AddMessage function rather than print to display the tool's messages in the results window.

GIS 5103 Module 9: Debugging and Error Handling

This week's module formally introduced us to the tools and tasks involved in debugging scripts.  While I'm sure we have all been forced to debug our scripts throughout the course, this module taught us how best to handle the inevitable errors that will plague us in our careers.  The debugging tool in PythonWin and try-except statements in particular will be useful in the future.

This assignment differed from previous ones in that we were given three scripts and tasked with debugging them.  The first two scripts contained several errors we needed to locate and correct for them to run properly.  For the final script, we were to insert a try-except statement that would trap the error in part A, print the error message to the interactive window, and continue to run part B of the script.  Below are screenshots of the interactive window messages resulting from my corrected scripts.

Script 1 interactive window after debugging.

Script 2 interactive window after debugging.

Script 3 interactive window after adding try-except statement.

GIS 5103 Module 8: Working with Rasters

This week saw us experiment with manipulating raster data using Python script.  As we did previously with vector data, we learned how to return information on raster files, such as spatial reference and cell size.  Because rasters are often large files, Python script will reference a raster object and create temporary rasters rather than create multiple raster files.  This is efficient, but it does require one to remember to save the final raster, otherwise nothing will remain after the script is complete.  

As in previous modules, our assignment required us to utilize many of the tools learned about in the exercise and expand them a bit further.  After creating a for loop to ensure the Spatial Analyst extension was available, our script needed to go through several steps.  First, it needed to use the RemapValue and Reclassify functions to collapse three separate classification of the landcover raster to one.  It then needed to create temporary rasters that would be used to create a final output between specific ranges of slope and aspect.  To do this, we needed to create four temporary rasters.  Two of the rasters were to define the slope range, with one including slope less than the maximum value while the other included slope above the minimum value.  The other two rasters would do the same for aspect.  These four rasters were to be combined, along with the reclassified landcover raster, using the "&" operator.  This would create a final raster with slope and elevation between the specified ranges along with the reclassified landcover data.  This final raster then needed to be saved in order to be permanent.  A screenshot of this raster can be seen above.  We were also required to insert print statements in our script to show in the interactive window the process of the script as it was run (see screenshot below).