3
source_srs = get_srs(shp_abs_path)
source_epsg = source_srs.GetAttrValue("PROJCS|GEOGCS|AUTHORITY", 1)
print source_epsg #prints 4269
target_epsg = 3857 # because these are headed to leafletjs
target_srs = osr.SpatialReference()
target_srs.ImportFromEPSG(target_epsg)
osr_transform = osr.CoordinateTransformation(source_srs, target_srs)

then as I enumerate through the shapefile I transform the bboxes on the fly

sf = shapefile.Reader(base_path)
shapeRecords = sf.shapeRecords()  #will store the geometry separately
self.shapefile_records[fn] = shapeRecords
for i, shape in enumerate(sf.shapes()):
    bbox = shape.bbox
    ## OGR/SRS TRANSFORMATION
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(bbox[0],bbox[1])
    point.Transform(osr_transform)
    print point.GetX(), point.GetY()

    point.AddPoint(bbox[2],bbox[3])
    point.Transform(osr_transform)
    print point.GetX(), point.GetY()

which prints coordinates like this...

-13613708.5846 4623834.1438
-13613668.296 4623869.60132

these put me out in the ocean I think, but they're supposed to be in Napa county. Any suggestions are warmly welcomed. I can use almost any python lib.

NEW

I've identified an issue and I'm not sure if it's with osgeo or the 4269 to 4326 transformation I've defined. Have a look at this output:

>>> from osgeo import ogr, osr
>>> source_srs  = osr.SpatialReference()
>>> source_srs.ImportFromEPSG(4269)
0
>>> target_srs = osr.SpatialReference()
>>> target_srs.ImportFromEPSG(4326)
0
>>> osr_transformation = osr.CoordinateTransformation(source_srs, target_srs)
>>> point = ogr.Geometry(ogr.wkbPoint)
>>> point.AddPoint(6449750.845, 2001628.41312)
>>> point.Transform(osr_transformation)
0
>>> print point.Transform(osr_transformation)
0
>>> point.GetX(), point.GetY()
(6449750.845, 2001628.4131199997)

Nothing is modified. I've triple confirmed the data is in EPSG:4269 and included a shortened bitly link to it below in the comments.

  • 1
    get_srs(shp_abs_path)come from ? source_srs is not an osr.SpatialReference() – gene Oct 19 '14 at 19:37
  • gene it returns an srs object from another routine - something along these lines bit.ly/1prDBWU but returns the srs object – Nona Yerbiz Oct 19 '14 at 19:44
  • It's not clear if you have done it, but it's a really good idea to use osr.UseExceptions() – Mike T Oct 19 '14 at 20:59
2

The problem with your solution is that source_srs is not a valid osr.SpatialReference(). If the result of source_epsg is 4269 then:

source_srs  = osr.SpatialReference()
source_srs.ImportFromEPSG(4269) 

gives a valid osr.SpatialReference()

NEW

If I understand correctly your question, you want to use Leaflet, and Leaflet expects coordinates and GeoJSON to be in EPSG:4326 (Reprojecting coordinates and geoJSON in leaflet)

If I project your results with osr

from osgeo import osr
wgs84 = osr.SpatialReference()
wgs84.ImportFromEPSG(4326)
merca  = osr.SpatialReference()
merca.ImportFromEPSG(3857)
transformation = osr.CoordinateTransformation(merca,wgs84)
print transformation.TransformPoint(-13613708.5846,4623834.1438)
-122.29402495095313, 38.313684256028715
transformation.TransformPoint(-13613668.296,4623869.60132)
-122.29366303230159, 38.313934175663768, 0.0

You can test the resulting GeoJSON:

 {"type":"FeatureCollection", "features": [
 { "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -122.29402495095313, 38.313684256028715 ] } },
 { "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -122.29366303230159, 38.313934175663768] } }
 ]}

In geosjon.io or GeJSONLint

enter image description here

The results are the same with EPSG:4269 (look at Silly Geographic Precision):

{"type":"FeatureCollection", "features": [  
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -122.29402495095315, 38.313684256946736 ] } },
{ "type": "Feature", "properties": { }, "geometry": { "type": "Point", "coordinates": [ -122.29366303230159, 38.313934176581782 ] } }
]}

With target_epsg = 3857:

enter image description here

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.