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.

  • geopandas is your friend for vector files. Rasterio for raster. – RutgerH Nov 30 '18 at 9:23

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.

enter image description here

Learning Geospatial Analysis with Python of Joel Lawhead

enter image description here

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)

  • What is this Python stuff of that you speak ;) – Nathan W Jun 22 '13 at 8:03
  • Fiona seems to be throwing a DLL error on Windows. – multigoodverse Apr 27 '16 at 14:42
  • How did you install Fiona ? no problem for me – gene Apr 27 '16 at 14:56

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

# 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'

# 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 = layer.GetNextFeature()

# close the data source and text file
  • 4
    .Destroy is an awesome method name :p – Jason Jun 21 '13 at 17:04

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.

  • 1
    please include enough of the essence of the solution can be grasped and understood before going to visit and read the page. In time that page will probably not be at that address rendering this answer not very useful. :) – matt wilkie Jun 21 '13 at 19:23

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

  • Any chance you could post some sample codes using shapely? – sherpas Jun 21 '13 at 16:33
  • 2
  • 5
    Link only answers aren't helpful in the long run, as they inevitably become broken. Please include enough information about the destination that a) it's new home can be rediscovered, and b) the essence of the solution can be grasped and understood before going to visit and read the page. – matt wilkie Jun 21 '13 at 19:20

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.