I do aviation maps for a living and basically most of our maps have a raster layer containing elevation. One of the things that I wanted to do recently is to check what is the maximum elevation in the view extent of the layout. There is a bit of a discussion in the twitterverse among #GIS people about what things you should know about as a GIS professional, one that comes high is the ability to code. I have taken a few lessons in python enough to be able to script a few things.

You should really learn python is you are into GIS, SQL & Javascript also!

MHLM

MHLM VOR Z RWY 22

Since our maps are done using ArcGIS, I decided to make a python script using the arcpy module to get the result I wanted. I was going to try initially to use model builder but for some reason I just went ahead and coded what I needed.

So what I needed was basically:

1. Get the active dataframe extent

2. Get the raster layer for clipping

3. Clip raster to active dataframe extent

4. Get Max and Min values for clipped raster

Caveats: dataframe and raster must be in same projection

[code language=”python”]
”’
This script requires data frame and raster in the same coordinate system
”’

# import python module
import arcpy

#get current mxd

mxd = arcpy.mapping.MapDocument(‘CURRENT’)

#get active dataframe

adf = mxd.activeDataFrame

# get extent coordinates
xmax = adf.extent.XMax
xmin = adf.extent.XMin
ymax = adf.extent.YMax
ymin = adf.extent.YMin

extent = ‘%s %s %s %s’ % (xmin,ymin,xmax,ymax)

#get raster layer

for lyrs in arcpy.mapping.ListLayers(mxd,"SRTM_FIR_MHTG_FEET_tinted",adf):
clip_raster = lyrs.name

#clip raster to data frame extent creating a in memory raster layer

arcpy.Clip_management(clip_raster,extent, "in_memory/elev_raster", "#", "#", "NONE")

#get elev_raster maximum value

elevSTDResultMax = arcpy.GetRasterProperties_management("in_memory/elev_raster", "MAXIMUM")
elevSTDResultMin = arcpy.GetRasterProperties_management("in_memory/elev_raster", "MINIMUM")
elevSTDMax = elevSTDResultMax.getOutput(0)
elevSTDMin = elevSTDResultMin.getOutput(0)
print elevSTDMin,elevSTDMax
[/code]

Since I used a memory layer to store the raster once we remove it and close the session it will no longer exist, you could write the result to disk but I have no need for that in my situation. So now just open the Python window and load the py file and voila …

arcpy clip result

arcpy clip result

Min: 23 ft | Max: 5653 ft

Now this is a basic script that I intend to build a few other things on top of it but you got to start somewhere.


0 Comments

Leave a Reply

Avatar placeholder

Your email address will not be published.

%d bloggers like this: