[GIS] Python Script examples for geoprocessing shapefiles without using arcpy

geoprocessingpython

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.

Best Answer

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 !):

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
(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)

Related Question