Proj4j is reprojecting EPSG:32614 differently than OGR

[lossyrob]

If I reproject this LatLng polygon:

POLYGON ((-96.952497 37.97709252761919, -96.952497 38.143939, -96.78565052761918 38.143939, -96.78565052761918 37.97709252761919, -96.952497 37.97709252761919))

from LatLng to CRS.fromName("EPSG:32614"), I get

POLYGON ((679830.5797888724 4205251.08962126, 679830.5797888724 4224102.15528983, 694044.9218645435 4224102.15528983, 694044.9218645435 4205251.08962126, 679830.5797888724 4205251.08962126))

If I reproject it with OGR, I get

POLYGON ((679830.5797888713 4205251.08962126, 679422.5041216067 4223766.106720305, 694044.9218645423 4224102.155289831, 694486.2944102848 4205586.653099069, 679830.5797888713 4205251.08962126))

These polygons are different. Why?

[dwins]

I'm not finding significant differences on my build. I wrote a small test harness to use the OGR python bindings and geotrellis' build of proj4j to reproject the sample input from this issue and print the absolute differences for each coordinate. The result I got was (1.0477378964424133E-9,0.0),(1.1641532182693481E-9,-9.313225746154785E-10),(1.1641532182693481E-9,-9.313225746154785E-10),(1.0477378964424133E-9,0.0),(1.0477378964424133E-9,0.0) ; much smaller than the reported discrepancy. I'm not sure what might have changed but the code and the epsg definitions seem pretty much identical between proj.4 and proj4j so I think we can put this one down to differences between how Python and Java are parsing floating point values from strings.

For the record, my test code looks like this:

import geotrellis.vector._                                                                                                                                                                                         
import geotrellis.vector.io.wkt._                                                                                                                                                                                  
import geotrellis.vector.reproject._                                                                                                                                                                               
import geotrellis.proj4._                                                                                                                                                                                          
import scala.sys.process._                                                                                                                                                                                         

object Main {                                                                                                                                                                                                      
  implicit class withAsStream(val s: String) extends AnyVal {                                                                                                                                                      
    def asStream: java.io.InputStream = new java.io.ByteArrayInputStream(s.getBytes("utf-8"))                                                                                                                      
  }                                                                                                                                                                                                                

  def main(args: Array[String]) {                                                                                                                                                                                  
    val pythonScript = "./check.py 4326 32614"                                                                                                                                                                     
    val p = """POLYGON ((-96.952497 37.97709252761919, -96.952497 38.143939, -96.78565052761918 38.143939, -96.78565052761918 37.97709252761919, -96.952497 37.97709252761919))"""                                 
    val p2 = p.parseWKT.asInstanceOf[Polygon].reproject(CRS.fromName("epsg:4326"), CRS.fromName("epsg:32614"))                                                                                                     
    val p3 = (pythonScript #< p.asStream !!).parseWKT.asInstanceOf[Polygon]                                                                                                                                        
    println((p2.vertices, p3.vertices).zipped.map { case (a, b) => (a.x - b.x, a.y - b.y) }.mkString(","))                                                                                                         
  }                                                                                                                                                                                                                
}

Accompanied by a Python script:

#!/usr/bin/env python                                                                                                                                                                                              
from osgeo import ogr                                                                                                                                                                                              
from osgeo import osr                                                                                                                                                                                              
import sys                                                                                                                                                                                                         

source = osr.SpatialReference()                                                                                                                                                                                    
source.ImportFromEPSG(int(sys.argv[1]))                                                                                                                                                                            

target = osr.SpatialReference()                                                                                                                                                                                    
target.ImportFromEPSG(int(sys.argv[2]))                                                                                                                                                                            

transform = osr.CoordinateTransformation(source, target)                                                                                                                                                           

wkt_text = sys.stdin.read()                                                                                                                                                                                        
poly = ogr.CreateGeometryFromWkt(wkt_text)                                                                                                                                                                         
poly.Transform(transform)                                                                                                                                                                                          

print poly.ExportToWkt()

[lossyrob]

Looks like this issue is a no-go, closing.