Python Script examples for geoprocessing shapefiles without using arcpy

Question

I would like to use a Python script that is not based on arcpy to do things like query a shapefile by attributes, create new layer from selection, and calculate areas of a polygon and convert polygons to points.

Anyone have any code examples of using other Python modules or libraries? I am able to do this easily using arcpy but i wanted to explore other options.

Answer 1

That's strange, as if people suddenly discovered the power of Python (without ArcPy which is just one Python module among others), see for example the question Visualize shapefile in Python:

• geospatial processing in Python has a very long history, much older than Arcpy (or arcgisscripting) -> no "mimic" the capabilities of ArcPy here, as Paul says, most were already there before ArcPy.
• the reference for the Python modules is the Python Package Index (Pypi) and there is a dedicated section: Topic :: Scientific/Engineering :: GIS
• you can do anything with these modules and it is often easier and faster than ArcPy because it is pure Python (no cursors...).
• Shapely is one of these modules for processing geospatial geometries -> calculate areas of a polygon and convert polygons to points..
• if you want to process vector layers, there is osgeo/ogr, Fiona or Pyshp (and others, less used) -> query a shapefile by attributes, create new layer from selection, calculate areas of a polygon, convert polygons to points
• for processing rasters, the standard is osgeo/gdal
• for spatial analysis, there is Pysal
• for 3D, you can use other Scientific modules like numpy or scipy (3D algorithms, grids, but also statistics, geostatistics, 2D or 3D)
• And I don't talk about mapnik, matplotlib/basemap,Geodjango and ...

You can combine all (Pysal with shapely, ...) and mix them with the other Scientific modules.

Thus for Python Script examples, search for Pyshp Fiona, ogr, gdal or shapely in gis.stackexchange or the internet (many examples, not only in English).)
One of them in French (the scripts and the figures are universal !):

• Python: Using vector and raster layers in a geological perspective, without GIS software
  an other in English:
• GIS with Python, Shapely, and Fiona
  and in Spanish
• Determination of areas of irregular polygons using the coordinates of the vertices
  in gis.stackexchange
• Elevation profile 10 km each side of a line
• Updating Attributes using Pyshp
• How to create a 3D shapefile from a raster?
• Python Script for getting elevation difference between two points
• etc

The script presented by Aaron can be written more simply with Fiona that uses only Python dictionaries:

import fiona
with fiona.open('sites.shp', 'r') as input:
    with open('hw1a.txt', 'w') as output:
       for pt in input:
           id = pt['properties']['id']
           cover = pt['properties']['cover']
           x = str(point['geometry']['coordinates'][0])
           y = str(point['geometry']['coordinates'][21])
           output.write(id + ' ' + x + ' ' + y+ ' ' + cover + '\n')

and if you use shapely in addition:

from shapely.geometry import shape
with fiona.open('sites.shp', 'r') as input:
    with open('hw1a.txt', 'w') as output:
       for pt in input:
           id = pt['properties']['id']
           cover = pt['properties']['cover']
           x = str(shape(pt['geometry']).x)
           y = str(shape(pt['geometry']).y)
           output.write(id + ' ' + x + ' ' + y+ ' ' + cover + '\n')

There are also two books:

Python Geospatial Development of Eric Westra.

Learning Geospatial Analysis with Python of Joel Lawhead

_{(source: cloudfront.net)}

Python is also used as a scripting language in other GIS applications like QGIS (Quantum GIS), GRASS GIS, gvSIG or OpenJump or 3D modelers like Paraview (and Blender also !). And you can use the majority of the geospatial modules in all these application (see Visualising QGIS data with Blender)

Answer 2

I highly recommend the USU site Geoprocessing with Python using Open Source GIS to get you started. They primarily use the GDAL/OGR library throughout the exercises. Installing GDAL/OGR can be a bit of a challenge, so this blog entry may be helpful for you: Installing GDAL (and OGR) for Python on Windows. Also check out Alternatives to using Arcpy on GIS.SE.

The following opensource geoprocessing script example (from the USU site) is used to extract attribute data and write it to a text file:

# import modules
import ogr, os, sys

# set the working directory
os.chdir('f:/data/classes/python/data')

# open the output text file for writing
file = open('hw1a.txt', 'w')

# get the shapefile driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# open the data source
datasource = driver.Open('sites.shp', 0)
if datasource is None:
  print 'Could not open file'
  sys.exit(1)

# get the data layer
layer = datasource.GetLayer()

# loop through the features in the layer
feature = layer.GetNextFeature()
while feature:

  # get the attributes
  id = feature.GetFieldAsString('id')
  cover = feature.GetFieldAsString('cover')

  # get the x,y coordinates for the point
  geom = feature.GetGeometryRef()
  x = str(geom.GetX())
  y = str(geom.GetY())

  # write info out to the text file
  file.write(id + ' ' + x + ' ' + y + ' ' + cover + '\n')

  # destroy the feature and get a new one
  feature.Destroy()
  feature = layer.GetNextFeature()

# close the data source and text file
datasource.Destroy()
file.close()

Answer 3

You might be interested in GDAL/OGR.

GDAL is used for processing rasters while OGR is used for vectors. Both are open source libraries.

If you are looking to remove some dependency on ArcPy, you can mimic some capabilities by reading the information to an array and running your own calculations in pure Python.

I recently did this with selecting points in a polygon, as seen here. It utilizes the ray casting algorithm to determine if a point lies within a polygon, given the coordinates of the polygon's vertices.

Answer 4

I've never used this personally, but others in the office like to use shapely: https://pypi.python.org/pypi/Shapely