Different definition of EPSG:3857 in GeoTools and QGIS

Question

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?

Answer 1

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"]]'

Answer 2

CRS is not part of the original ESRI shapefile format. QGIS uses two rivaling extension files .qpj and .prj to store CRS information. Both are text files containing Well Known Text (WKT).

QPJ is QGIS specific, as far as I know it is not read by other software.

PRJ was initially defined by ESRI and exists in various flavours. I assume QGIS tries to create maximal compatibility. PRJ files from recent QGIS versions do not contain any EPSG reference numbers. I think GeoTools reads them as custom projection and doesn't link it to any predefined CRS.

See also Is there a standard for the specification of .prj files?