projection-latlng-wm-differences.scala

/** GDAL will translate a GeoTIFF with these boundies in EPSG:4326 -
  *    Upper Left  ( -85.4969788,  35.4371767) ( 85d29'49.12"W, 35d26'13.84"N)
  *    Lower Left  ( -85.4969788,  35.3464955) ( 85d29'49.12"W, 35d20'47.38"N)
  *    Upper Right ( -85.2701788,  35.4371767) ( 85d16'12.64"W, 35d26'13.84"N)
  *    Lower Right ( -85.2701788,  35.3464955) ( 85d16'12.64"W, 35d20'47.38"N)
  *    Center      ( -85.3835788,  35.3918361) ( 85d23' 0.88"W, 35d23'30.61"N)
  * 
  * to these boundry points in EPSG:3857 -
  *    Upper Left  (-9517480.140, 4223451.561) ( 85d29'49.12"W, 35d26'13.84"N)
  *    Lower Left  (-9517480.140, 4211063.427) ( 85d29'49.12"W, 35d20'47.24"N)
  *    Upper Right (-9492230.441, 4223451.561) ( 85d16'12.56"W, 35d26'13.84"N)
  *    Lower Right (-9492230.441, 4211063.427) ( 85d16'12.56"W, 35d20'47.24"N)
  *    Center      (-9504855.291, 4217257.494) ( 85d23' 0.84"W, 35d23'30.59"N)
  * 
  * (GDAL command: gdalwarp -t_srs EPSG:3857 latlng.tiff webmerc.tiff)
  * 
  * ogr2ogr will take this GeoJSON:
  *   {
  *    "type": "Feature",
  *    "properties": {},
  *    "geometry": {
  *      "type": "Polygon",
  *      "coordinates": [
  *         [
  *           [
  *             35.43717666390691, -85.26763916015625
  *           ],
  *           [
  *             35.43717666390691, -85.49697875976562
  *           ],
  *           [
  *             35.34649548039281, -85.49697875976562
  *           ],
  *           [
  *             35.34649548039281, -85.26763916015625
  *           ],
  *           [
  *             35.43717666390691, -85.26763916015625
  *           ]
  *         ]
  *       ]
  *     }
  *   }
  * 
  * and transform it from EPSG:4326 to EPSG:3857 to this:
  *   {
  *     "type": "FeatureCollection",
  *     "crs": {
  *       "type": "name",
  *       "properties": {
  *         "name": "urn:ogc:def:crs:EPSG::3857"
  *       }
  *     },
  *     "features": [
  *       {
  *         "type": "Feature",
  *         "properties": {},
  *         "geometry": {
  *           "type": "Polygon",
  *           "coordinates": [
  *             [
  *               [
  *                 3944848.4613773944,
  *                 -20323175.971143067
  *               ],
  *               [
  *                 3944848.4613773944,
  *                 -20640357.197044015
  *               ],
  *               [
  *                 3934753.8782040733,
  *                 -20640357.197044015
  *               ],
  *               [
  *                 3934753.8782040733,
  *                 -20323175.971143067
  *               ],
  *               [
  *                 3944848.4613773944,
  *                 -20323175.971143067
  *               ]
  *             ]
  *           ]
  *         }
  *       }
  *     ]
  *   }
  * 
  * (ogr2ogr command: ogr2ogr -s_srs EPSG:4326 -t_srs EPSG:3857 polygon_wm.json polygon.json -f GeoJSON)
  * 
  * GeoTools will transform this Polygon in EPSG:4326 -
  *    POLYGON ((35.43717666390691 -85.26763916015625, 
  *              35.43717666390691 -85.49697875976562, 
  *              35.34649548039281 -85.49697875976562, 
  *              35.34649548039281 -85.26763916015625, 
  *              35.43717666390691 -85.26763916015625))
  * 
  * to this Polygon in EPSG:3857 -
  *    POLYGON ((-9491950.172453186 4223451.560869129, 
  *              -9517480.139900435 4223451.560869129, 
  *              -9517480.139900435 4211068.7622869285, 
  *              -9491950.172453186 4211068.7622869285, 
  *              -9491950.172453186 4223451.560869129))
  * 
  * Proj4J will transform the same Polygon in EPSG:4326 to this Polygon in EPSG:3857 -
  *    POLYGON ((3944848.4613773944 -20323175.97114305, 
  *              3944848.4613773944 -20640357.197043996, 
  *              3934753.8782040733 -20640357.197043996, 
  *              3934753.8782040733 -20323175.97114305, 
  *              3944848.4613773944 -20323175.97114305))
  * 
  * See code below for how this was generated.
  * 
  */

import com.vividsolutions.jts.geom._

val factory = new GeometryFactory()

val lr =
  factory.createLinearRing(
    Array(
      new Coordinate(35.43717666390691, -85.26763916015625),
      new Coordinate(35.43717666390691, -85.49697875976562),
      new Coordinate(35.34649548039281, -85.49697875976562),
      new Coordinate(35.34649548039281, -85.26763916015625),
      new Coordinate(35.43717666390691, -85.26763916015625)
    )
  )

val polygon = factory.createPolygon(lr, Array[LinearRing]())

def reprojectProj4J(polygon: Polygon): Polygon = {
  import org.osgeo.proj4j._

  val crsFactory = new CRSFactory()
  val latLng = crsFactory.createFromName("EPSG:4326")
  val webMercator = crsFactory.createFromName("EPSG:3857")
  
  val transformFactory = new CoordinateTransformFactory()
  val t = transformFactory.createTransform(latLng, webMercator)

  val coords: Array[Coordinate] = polygon.getCoordinates
  val projectedCoords =
    coords.map { coord =>
      val srcP = new ProjCoordinate(coord.x, coord.y)
      val destP = new ProjCoordinate()

      t.transform(srcP, destP)
      new Coordinate(destP.x, destP.y)
    }
  val projectedLinearRing = factory.createLinearRing(projectedCoords)
  factory.createPolygon(projectedLinearRing, Array[LinearRing]())
}

def reprojectGeoTools(polygon: Polygon): Polygon = {
  import org.geotools.geometry.jts.JTS
  import org.geotools.referencing.CRS

  val latLng = CRS.decode("EPSG:4326")
  val webMercator = CRS.decode("EPSG:3857")

  val transform = CRS.findMathTransform(latLng, webMercator, true)
  val transformed = JTS.transform(polygon, transform).asInstanceOf[Polygon]
  transformed
}

val proj4Poly = reprojectProj4J(polygon)
val geotoolsPoly = reprojectGeoTools(polygon)

println(polygon)
println(proj4Poly)
println(geotoolsPoly)
proj4Poly should be (geotoolsPoly)