Category Archives: Python

Book Review: Python Scripting for ArcGIS by Paul A. Zandbergen

Posted on by

Title: Python Scripting for ArcGIS
Author: Paul A. Zanbergen
Publisher: ESRI Press
Year: 2013
Aimed at: Python/ArcPy – beginners, ArcGIS – knowledgeable
Purchased from: www.bookdepository.com

Python Scripting for ArcGIS

This book is a fantastic stepping stone for beginners into the enchanted world of ArcPy. ArcPy is a Python site package that provides access to the extensive set of geoprocessing tools available in ArcGIS. Besides enabling programmatic geospatial analysis ArcPy modules also facilitate data management, data conversion and map document management.

I think a quote from the Preface pages of this book aptly sums up what the book is all about.

“a little bit of code goes a long way.”

As an introductory text your eyes will be opened to how small snippets of code can run geoprocessing tools that can form the basis for extensive geospatial analysis. You won’t find in-depth spatial analysis or data management techniques but you will find an easy to read, easy to follow informative text book that provides the theory behind using Python/ArcPy and will act as a reference to the capabilities of ArcPy.

Before purchasing this book I read a number of reviews. While an overwhelming majority applauded the book there where a few who complained about the basic introduction to Python provided. Even though there is a chapter dedicated to creating Python functions and classes one review that sticks out in my mind wanted in-depth object orientated programming for GIS Python which to me is miles beyond the scope of this book. The author does a great job of providing a primer to the Python language but this is not what this book is about. There are a myriad of Python text books for beginners and also online tutorials out there and I would certainly recommend making use of these and getting comfortable with the general syntax, data structures and data types before diving head first into using Python for geospatial activities.

I bought this book because I wanted a foundation for ArcPy that I could build upon. While progressing through the text I was constantly looking to the ArcGIS Resources pages for more information about geoprocessing tools encountered and the syntax required to implement them programmatically. I would recommend using this book in tandem with the Resource pages for the ultimate beginner experience. The book is extremely informative for a beginner’s text but it will be your genuine interest in the material that will take you well beyond what’s on offer here.

The book and topics are well designed with each chapter building upon the previous. The first part introduces the Python language, development environments (PythonWIn and the Interactive Python WIndow in ArcMap), and the basics of geoprocessing. Part two is where you begin your ArcPy experience, writing scripts and learning about ArcPy modules and their capabilities. Part three introduces some specialized tasks such as automating ArcMap workflows through map scripting and error handling is also discussed. Part four provides an introduction to creating your own custom tool.

Some of the more interesting materials I found covered in this book were; working with the mapping module for automating map document tasks, accessing and manipulating data with cursors and the data access module, working with geometries and rasters, and creating custom tools. These will provide the springboard for you to dive into more advanced scripting.

Overall Verdict: The book was a great investment (c. €60). It would be hard to find a better way to introduce yourself to ArcPy. It won’t teach you everything you need to know to build applicable scripts but provides an invaluable foundation. Highly recommended for beginners.

Generate a Projection File (.prj) using Python

Posted on by

This blog has moved to a new home at Final Draft Mapping.

Generate a Projection File (.prj) using Python

A .prj file contains the coordinate system information for a dataset and is required for ‘on the fly’ projection by GIS software. The file itself is a text file containing information in Well Known Text (WKT) format. The code snippet here shows you how to generate a .prj file, using Python, for your data if the .prj file is missing. See the post on CSV to Shapefile with pyshp for an example of using the function in a workflow. You will need to know what coordinate system the data is in and the EPSG code. For example, to find the EPSG code for WGS84 use a search engine and search for EPSG WGS84. The first site returned is usually http://spatialreference.org/. The EPSG code for WGS84 is 4326.

# function to generate .prj file information using spatialreference.org
def getWKT_PRJ (epsg_code):
     import urllib
     # access projection information
     wkt = urllib.urlopen("http://spatialreference.org/ref/epsg/{0}/prettywkt/".format(epsg_code))
     # remove spaces between charachters
     remove_spaces = wkt.read().replace(" ","")
     # place all the text on one line
     output = remove_spaces.replace("\n", "")
     return output

# create the .prj file
prj = open("filename.prj", "w")
# call the function and supply the epsg code
epsg = getWKT_PRJ("4326")
prj.write(epsg)
prj.close()

The filename should match the name of your data files. For examples, the .prj file for cities.shp or cities.tab should be cities.prj.

If you open the .prj file in a text editor you will see the text below.

GEOGCS["WGS84",DATUM["WGS_1984",SPHEROID["WGS84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

GEOGCS at the beginning of the WKT indicates that EPSG: 4326 is a geographic coordinate reference system (it would be PROJCS if it was a projected coordinate system). The information within the square brackets after DATUM provide the information for the parameters of the datum. We see the name of the datum, ‘WGS_1984’, and the SPHEROID information with a semimajor axis of 6,378,137 m with an inverse-flattening ratio of 298.257223563. The PRIMEM data tells us that it uses Greenwich as the prime meridian (where longitude is zero). UNIT specifies the measurement unit of the coordinate system, here it is ‘degree’, and the 0.0174532925199433 value is required to convert from radians into the stated units.

Open your data in a GIS and make sure that it is displaying correctly in relation to other features.

Decimal-Degrees (DD) to Degrees-Minutes-Seconds (DMS) with Python

Posted on by
3

This blog post has moved to a new home at Final Draft Mapping

https://learn.finaldraftmapping.com/convert-decimal-degrees-dd-to-degrees-minutes-seconds-dms-with-python/

I am currently transitioning from the comfort of hiding behind buttons and converting to using programming for GIS and geospatial & geostatistical analysis. The code below is an implementation of converting decimal degrees to degrees-minutes-seconds format. The code is heavily commented and hopefully easy to understand.

import math

def dd2dms(longitude, latitude):

    # math.modf() splits whole number and decimal into tuple
    # eg 53.3478 becomes (0.3478, 53)
    split_degx = math.modf(longitude)
    
    # the whole number [index 1] is the degrees
    degrees_x = int(split_degx[1])

    # multiply the decimal part by 60: 0.3478 * 60 = 20.868
    # split the whole number part of the total as the minutes: 20
    # abs() absoulte value - no negative
    minutes_x = abs(int(math.modf(split_degx[0] * 60)[1]))

    # multiply the decimal part of the split above by 60 to get the seconds
    # 0.868 x 60 = 52.08, round excess decimal places to 2 places
    # abs() absoulte value - no negative
    seconds_x = abs(round(math.modf(split_degx[0] * 60)[0] * 60,2))

    # repeat for latitude
    split_degy = math.modf(latitude)
    degrees_y = int(split_degy[1])
    minutes_y = abs(int(math.modf(split_degy[0] * 60)[1]))
    seconds_y = abs(round(math.modf(split_degy[0] * 60)[0] * 60,2))

    # account for E/W & N/S
    if degrees_x < 0:
        EorW = "W"
    else:
        EorW = "E"

    if degrees_y < 0:
        NorS = "S"
    else:
        NorS = "N"

    # abs() remove negative from degrees, was only needed for if-else above
    print "\t" + str(abs(degrees_x)) + u"\u00b0 " + str(minutes_x) + "' " + str(seconds_x) + "\" " + EorW
    print "\t" + str(abs(degrees_y)) + u"\u00b0 " + str(minutes_y) + "' " + str(seconds_y) + "\" " + NorS

# some coords of cities
coords = [["Dublin", -6.2597, 53.3478],["Paris", 2.3508, 48.8567],["Sydney", 151.2094, -33.8650]]

# test dd2dms() 
for city,x,y in coords:
    print city + ":"
    dd2dms(x, y)

The output from running the above script is…

Dublin:
    6° 15' 34.92" W
    53° 20' 52.08" N
Paris:
    2° 21' 2.88" E
    48° 51' 24.12" N
Sydney:
    151° 12' 33.84" E
    33° 51' 54.0" S

I am happy that the output of this script is the same as online converters.