

In [1]:

    
name = '2016-06-10-arcgis-intro'
title = 'Introduction to ArcGIS and its Python interface'
tags = 'gis, maps, basics'
author = 'Melanie Froude'



In [2]:

    
from nb_tools import connect_notebook_to_post
from IPython.core.display import HTML

html = connect_notebook_to_post(name, title, tags, author)

Today Melanie lead the meeting with a session on the ArcGIS software and how we can use Python to automatise the geospatial data processing. The slides are available below.

We started with a brief introduction to the types of data and analysis you can do in ArcGIS. Then Melanie demonstrated how to produce a 3D terrain model using the ArcScene toolbox.

Presentation



In [3]:

    
# embed pdf into an automatically resized window (requires imagemagick)
w_h_str = !identify -format "%w %h" ../pdfs/arcgis-intro.pdf[0]
HTML('<iframe src=../pdfs/arcgis-intro.pdf width={0[0]} height={0[1]}></iframe>'.format([int(i)*0.8 for i in w_h_str[0].split()]))









    Out[3]:

We all agreed that ArcGIS has a lot to offer to geoscientists. But what makes this software even more appealing is that you can work in a command-line interface using Python (ArcPy module).

So we looked at how to run processes using the Python window command-by-command and how you might integrate ArcGIS processes within a longer script. This was exemplified by Melanie's script that she used to analyse vegetation regrowth after a volcanic eruption.

The script takes two vegetation photos in GeoTIFF format retrieved by Landsat as input and calculates the Normalised Difference Vegetation Index (NDVI) for each of them. We can then compare the output to see how vegetation has changed over the time period.

Standard ArcGIS uses Python 2.7 (Python 3 is available in ArcGIS Pro)

The commands below require ArcGIS installed, and hence are not in executable cells in this notebook. </div>

ArcPy script example: NDVI of the two geotiff images

Import modules

import arcpy, string, arcpy.sa
from arcpy import env

Check out extension and set overwrite outputs

arcpy.CheckOutExtension("spatial")
arcpy.env.overwriteOutput = True

Stop outputs being added to the map

arcpy.env.addOutputsToMap = "FALSE"

Set workspace and declare variations

env.workspace = ("/path/to/demo/demo1")
print(arcpy.env.workspace)

Load the data

rasterb3 = arcpy.Raster("p046r28_5t900922_nn3.tif")
rasterb4 = arcpy.Raster("p046r28_5t900922_nn4.tif")

Describe variables

desc = arcpy.Describe(rasterb4)
print(desc.dataType)
print(desc.meanCellHeight)

Calculate the NDVI

Num = arcpy.sa.Float(rasterb4-rasterb3)
Denom = arcpy.sa.Float(rasterb4 + rasterb3)
NDVI1990 = arcpy.sa.Divide(Num, Denom)

Save the result as another .tif image

NDVI1990.save("/path/to/demo/demo1/NDVI1990.tif")

Do the same calculation for the images from a later year

rasterb3a = arcpy.Raster("L71046028_02820050721_B30.TIF")
rasterb4a = arcpy.Raster("L71046028_02820050721_B40.TIF")

Num = arcpy.sa.Float(rasterb4a-rasterb3a)
Denom = arcpy.sa.Float(rasterb4a + rasterb3a)
NDVI2005 = arcpy.sa.Divide(Num, Denom)

And after saving the second result, calculate the NDVI difference

NDVI2005.save("/path/to/demo/demo1/NDVI2005.tif")

NDVIdiff = NDVI2005 - NDVI1990
NDVIdiff.save("/path/to/demo/demo1/NDVIdiff.tif")

The result is shown in the slide 5.



In [4]:

    
HTML(html)









    Out[4]:





    
     This post was written as an IPython (Jupyter) notebook. You can view or download it using
    nbviewer.