python - Shapefile projection for Matplotlib Basemap

Preface: First time working with GIS data.

I'm using python's Basemap library for overlaying data on maps.

The tutorial I'm using includes a shapefile for London (london_wards.shp).

I use the bounds method to get the minimum bounding rectangle.

s="/data/london_wards.shp"
shp = fiona.open(s)
shp.bounds

(-0.5103750689005356, 51.28676016315085, 0.3340155643740321, 51.691874116909894)

This is the latitude, longitude (51.5072° N, 0.1275° W) for London.

I pass these bounds to a basemap instance.

---

Next, I've done the same for shapefiles obtained from SF OpenData.

However, the bounds are apparently not latitude, longitude.

s="/data/sfpd_districts.shp"
shp = fiona.open(s)
shp.bounds

(5979385.3163340045, 2085856.3243659574, 6024751.466714004, 2131290.9014959573)

So, I can't pass them to basemap.

Several shapefiles from SF OpenData have similar units returned from shp.bounds.

What are these values?

How are they converted to latitude, longitude for further use?

Answer

Use directly the GDAL module rather than using it through (Geo)Django:

from osgeo import osr
ori =  osr.SpatialReference()
desti =  osr.SpatialReference()
ori.ImportFromProj4("+init=EPSG:2227")
desti.ImportFromProj4("+init=EPSG:4326")
transformation = osr.CoordinateTransformation(ori,desti)
transformation1.TransformPoint(5979385.3163340045, 2085856.3243659574)
(-122.5129551683127, 37.706167841402952, 0.0)

or the pyproj module (with preserve_units=True because pyproj assumes that your coordinates are in meters and this is not the case for EPSG:2227):

from pyproj import Proj, transform
ori = Proj(init='epsg:2227',preserve_units=True) 
dest= Proj(init='EPSG:4326',preserve_units=True)
transform(ori, dest,x,y)
(-122.51295516831273, 37.706167841402959)

Combine with Fiona:

from pyproj import Proj, transform
import fiona

from fiona.crs import from_epsg
with fiona.open("/data/sfpd_districts.shp") as shp
    ori = Proj(shp.crs,preserve_units=True ),
    dest= Proj(init='EPSG:4326',preserve_units=True)
    with fiona.open('/data/sfpd_districtsWGS84.shp', 'w', 'ESRI Shapefile', shp.schema.copy(), crs=from_epsg(4326) as output:
        for point in shp:
            x,y =  point['geometry']['coordinates']
            point['geometry']['coordinates'] = transform(ori, dest,x,y)
            output.write(point)