[GIS] Different definition of EPSG:3857 in GeoTools and QGIS

coordinate systemgeotoolsqgis

I want to do some checking if my point lies in any of the polygons defined in shapefile. So, I thought that I will determine CRS from shapefile (created in QGIS), and use it for transformation (where source CRS is EPSG:4326). My shapefile has geometry in EPSG:3857 (WGS-87 pseudo mercator). However, when I get CRS programatically from shapefile such way (lets assume that featureSource is non null FeatureSource type object):

CoordinateReferenceSystem crs = featureSource.getSchema().getCoordinateReferenceSystem();

This result:
PROJCS["WGS_84_Pseudo_Mercator",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984", 6378137.0, 298.257223563]],
PRIMEM["Greenwich", 0.0],
UNIT["degree", 0.017453292519943295],
AXIS["Longitude", EAST],
AXIS["Latitude", NORTH]],
PROJECTION["Mercator_2SP"],
PARAMETER["standard_parallel_1", 0.0],
PARAMETER["latitude_of_origin", 0.0],
PARAMETER["central_meridian", 0.0],
PARAMETER["false_easting", 0.0],
PARAMETER["false_northing", 0.0],
UNIT["m", 1.0],
AXIS["x", EAST],
AXIS["y", NORTH]]

And when I check WKT definition of this CRS, it indeed is WGS87 pseudo mercator, but it is different from definition which I get that way:

CoordinateReferenceSystem crs = CRS.decode( "EPSG:3857" );

The result:
PROJCS["WGS 84 / Pseudo-Mercator",
GEOGCS["WGS 84",
DATUM["World Geodetic System 1984",
SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
UNIT["degree", 0.017453292519943295],
AXIS["Geodetic latitude", NORTH],
AXIS["Geodetic longitude", EAST],
AUTHORITY["EPSG","4326"]],
PROJECTION["Popular Visualisation Pseudo Mercator", AUTHORITY["EPSG","1024"]],
PARAMETER["semi_minor", 6378137.0],
PARAMETER["latitude_of_origin", 0.0],
PARAMETER["central_meridian", 0.0],
PARAMETER["scale_factor", 1.0],
PARAMETER["false_easting", 0.0],
PARAMETER["false_northing", 0.0],
UNIT["m", 1.0],
AXIS["Easting", EAST],
AXIS["Northing", NORTH],
AUTHORITY["EPSG","3857"]]

I assume that QGIS has different definition of it. What is also weird that if I use definition from the shapefile, my transformed geometry has incorrect coordinates, and if I use the one from geotools, coordinates are computed correctly.

But still, the question is, is geotools able to correctly determine CRS which is used with given shapefile? I know that QGIS is different story, but do you know what these differences come from?

Best Answer

You are using two methods from GeoTools for saying that there are different definitions of EPSG:3857 in GeoTools and QGIS. However, GeoTools and QGIS are both OSGeo Projects. You should use a QGIS method for checking the complete WKT definition of this CRS in QGIS. In this way, in QGIS:

crs_3857_WKT = QgsCoordinateReferenceSystem()
crs_3857_WKT.createFromSrid(3857)
True
crs_3857_WKT.toWkt()

resulting in a complete definition of this CRS in WKT and it obtained in QGIS (where their extension parameters are expressed in proj4 format).

u'PROJCS["WGS 84 / Pseudo-Mercator",GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]],PROJECTION["Mercator_1SP"],PARAMETER["central_meridian",0],PARAMETER["scale_factor",1],PARAMETER["false_easting",0],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AXIS["X",EAST],AXIS["Y",NORTH],EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs"],AUTHORITY["EPSG","3857"]]'

Related Question