from Python Script examples for geoprocessing shapefiles without using arcpy
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
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)