GeoTools Shapefile – Read Shapefile Without EPSG in PRJ File

epsggeotoolsshapefile

I'd like to improve the accuracy of reading a shapefile (password is "helpme") which has no EPSG specified in prj file.When I read the shapefile with geotools 21.0 I get this result (reprojected to EPSG:3857):

POLYGON ((-9009004.79558 4823535.01838, -9008959.57149 4823594.56743, -9008967.19005 4823606.26696, -9008958.6486 4823614.74616, -9008945.37693 4823608.60905, -9008929.96268 4823600.85573, -9008922.79142 4823597.16586, -9008874.70649 4823557.51187, -9008845.20358 4823534.47221, -9008822.41676 4823517.27486, -9008811.75253 4823499.34209, -9008875.43856 4823415.15485, -9008914.5396 4823448.50908, -9008945.19703 4823471.34222, -9008997.07237 4823514.39293, -9009002.26199 4823521.31009, -9009005.80562 4823524.12776, -9009008.31574 4823530.5046, -9009004.79558 4823535.01838))

Correct EPSG would be 32067. When I know this I get a better transformation with

MathTransform transform = CRS.findMathTransform(this.sourceCRS, this.targetCRS, true);

and then the result in targetCRS=3857 is:

POLYGON ((-9008970.13992 4823838.31668, -9008924.91457 4823897.8674, -9008932.53335 4823909.56725, -9008923.99165 4823918.0467, -9008910.71961 4823911.90942, -9008895.30493 4823904.15588, -9008888.13346 4823900.46592, -9008840.04717 4823860.81082, -9008810.54342 4823837.77053, -9008787.75595 4823820.5727, -9008777.09141 4823802.63943, -9008840.77923 4823718.44982, -9008879.88137 4823751.80497, -9008910.53968 4823774.63874, -9008962.41649 4823817.69065, -9008967.60626 4823824.608, -9008971.14999 4823827.42575, -9008973.66019 4823833.80277, -9008970.13992 4823838.31668))

which is far better (checked with OSM background layer).

When I open the shapefile with QGIS I get a dialog asking which transformation I want to use and just using the first proposed results to a better import when checking it against OSM map.
QGIS import result with OSM background layer

So if QGIS manages, how could I do better with GeoTools? I'm thinking of an utility that searches all EPSG definitions and finds the nearest to the params in my prj file.

Best Answer

GeoTools is reading your projection file just fine (though it would be best to ask your data provider to fix the prj file to use a standard format rather than one they just thought up over lunch one day). The problem is that you haven't specified which transform you want to use.

Versions of GeoTools released since GEOT-6920 was applied allow you to specify the transform you would like to use (rather like QGIS does). Since there hasn't been a release since the patch was added you will need to switch to the nightly SNAPSHOT builds of 26.x (main) currently to get this code to work.

MathTransform transform = CRS.findMathTransform( this.sourceCRS, this.targetCRS, "EPSG:8609");

If you are in Mississippi, other wise look up the transform code at https://epsg.io/32067 or you can look at the results of:

CoordinateReferenceSystem blmNad = CRS.decode("EPSG:32067");
CoordinateReferenceSystem webMerc = CRS.decode("EPSG:3857");
Map<String, AbstractCoordinateOperation> grid_trans = CRS.getTransforms(blmNad, webMerc);
for (String code : grid_trans.keySet()) {
  System.out
      .println(code + " " + grid_trans.get(code).getAccuracy());
}
Related Question