Skip to content

Instantly share code, notes, and snippets.

@lossyrob
Created July 15, 2014 21:25
Show Gist options
  • Save lossyrob/d997140f09425800028d to your computer and use it in GitHub Desktop.
Save lossyrob/d997140f09425800028d to your computer and use it in GitHub Desktop.
/** 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment