0

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?

1

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)

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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