3

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?

  • I don't understand your issue because I suppose that 'getCoordinateReferenceSystem' and 'decode' are methods of GeoTools. I can't see them in QgsCoordinateReferenceSystem of QGIS. Where is the comparison? – xunilk Mar 11 '15 at 22:08
  • This part: featureSource.getSchema().getCoordinateReferenceSystem(); returns CRS saved in shp file, and it comes from QGIS. This definition is different from definition from GeoTools (gt-epsg-hsql) which I get that way: CRS.decode( "EPSG:3857" ); – Daniel Szymatowicz Mar 12 '15 at 7:00
  • what is the output of the two methods? – Ian Turton Mar 12 '15 at 9:06
  • I updated the question with result WKTs. Thank you for taking a look. – Daniel Szymatowicz Mar 12 '15 at 20:30
  • 1
    The official definition is: epsg-registry.org/export.htm?wkt=urn:ogc:def:crs:EPSG::3857 – nmtoken May 10 at 15:47
2

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

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?

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.