Take the 2-minute tour ×
Geographic Information Systems Stack Exchange is a question and answer site for cartographers, geographers and GIS professionals. It's 100% free, no registration required.

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.

share|improve this question
add comment

4 Answers

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)

share|improve this answer
    
What is this Python stuff of that you speak ;) –  Nathan W Jun 22 '13 at 8:03
add comment

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()
share|improve this answer
3  
.Destroy is an awesome method name :p –  Jason Jun 21 '13 at 17:04
1  
+1 for actually including an example –  BenjaminGolder Jun 21 '13 at 20:58
add comment

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.

share|improve this answer
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
add comment

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

share|improve this answer
    
Any chance you could post some sample codes using shapely? –  user19351 Jun 21 '13 at 16:33
2  
4  
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
add comment

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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