OpenSource GIS Packages

ArcGIS is an OK solution for general GIS functions if it is available or you can afford to purchase it. QGIS is an open source GIS application that is included with the OSGeo4W (OSGeo for Windows) release. At the time of writing this (Spring 2015) OSGeo4W installs with a version of Python 2.7 but it does not have the correct file structure for Python to see "site-packages" so QGIS is not available. This leaves us without an open source general GIS library for Python at this time.

There are a large number of packages available that include GIS functionality. A portion of these are viewable on the Python.org web site.

OSGeo has taken over management of the core open source GIS toolkits.

Available Tools

The OpenSource tools include:

• GDAL - Raster file IO
• OGR - Vector file IO
• Proj4 - Projecting spatial data

Another point is that only GDAL comes with examples for Python scripting. OGR and Proj4 have Python bindings but without documentation (although this could change in the future). Another approach with these packages is to run the executable programs they provide as "subprocesses" which is well documented and heavily tested.

GeoPandas

As of 2018, GeoPandas appeared to be the best open source solution for accessing GIS functions in Python.

GeoPandas is a large set of libraries that include spatial functions. Documentation for GeoPandas is available on the GeoPandas web site. Below is some code that I used to begin evaluation of GeoPands to include it in GSP 318.

Note: GeoPandas is not on the lab computers and Brent is working on getting on a VM. I have had little luck installing GeoPandas directly. I was able to install Anaconda (which is huge!) and then install GeoPandas from Anaconda.

# import the GeoPandas libraries
import gdal
import geopandas
import geopandas as gpd

# import numpy and scipy
import numpy as np
from scipy import ndimage
import matplotlib.pylab as pylab
import matplotlib.pyplot as pyplot

# load a couple of shapefiles
world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))
cities = gpd.read_file(gpd.datasets.get_path('naturalearth_cities'))

# print the first few features of the file
print(world.head())

# GeoDataFrames are organized as named columns and then indexed arrays
# so this gets the number of entries in the 'name' column.
NumFeatures=len(world.name)

# access the attributes in the file to print all names
Count=0
while (Count<NumFeatures):
	print(world.name[Count])
	Count+=1

# print the bounds of the entire world
print(world.bounds)

# remove Antarctica from the world shapefile and then project it
world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]
world = world.to_crs({'init': 'epsg:3395'}) # world.to_crs(epsg=3395) would also work

world.plot() # create a map

pyplot.show() # show the map

print(world.crs)

# plot by GDP per cap

world = world[(world.pop_est>0) & (world.name!="Antarctica")]

world['gdp_per_cap'] = world.gdp_md_est / world.pop_est

world.plot(column='gdp_per_cap')

pyplot.show()

Additional Resources (outdated)

GDAL/OGR in Python

GDAL Tutorials (includes Python examples)