Shapefile projection for Matplotlib Basemap

Question

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?

gene · Accepted Answer · 2015-02-18T17:16:45.333

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)

Shapefile projection for Matplotlib Basemap

1 Answers1

Linked