How to parse CRS from shapefile using geotools

Question

I need to retrieve CRS from a shapefile. I have tried to get it using getCoordinateReferenceSystem() which is provided by org.geotools. My code is as follows.

URL someURL = new URL("https://github.com/basilmohammed/CRS-Extraction-From-ShapeFile/raw/master/ShapeFiles/3857/3857-Point/3857Point.shp");
ShapefileDataStore SHPdataStore = new ShapefileDataStore(someURL);
ContentFeatureSource featureSource = SHPdataStore.getFeatureSource();
ContentFeatureCollection featureCollection = featureSource.getFeatures();
SimpleFeatureIterator iterator = featureCollection.features();
SimpleFeatureType schema = featureSource.getSchema();
CoordinateReferenceSystem sourceCRS = schema.getCoordinateReferenceSystem();
// the coordinates system I want to reproject the data to
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
//Define a MathTransform object
MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);

But always I am getting sourceCRS with null as follows.

Any advice?

Answer 1

I was going to say that it was because you were opening the .prj file instead of the .shp file but surprisingly it works (because you are creating the ShapefileDatastore directly). The code you provided works just fine for me, except that it appears that, that shapefile is actually in WGS84 (EPSG:4326) already so the .prj file is wrong.

My Code

import java.io.IOException;
import java.net.URL;

import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.store.ContentFeatureCollection;
import org.geotools.data.store.ContentFeatureSource;
import org.geotools.geometry.jts.JTS;
import org.geotools.referencing.CRS;
import org.locationtech.jts.geom.Geometry;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.geometry.MismatchedDimensionException;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.NoSuchAuthorityCodeException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;

public class FetchCRS {

  public static void main(String[] args) throws IOException, NoSuchAuthorityCodeException, FactoryException,
      MismatchedDimensionException, TransformException {
    URL someURL = new URL(
        "https://github.com/basilmohammed/CRS-Extraction-From-ShapeFile/raw/master/ShapeFiles/3857/3857-Point/3857Point.shp");
    ShapefileDataStore SHPdataStore = new ShapefileDataStore(someURL);
    System.out.println(SHPdataStore);
    ContentFeatureSource featureSource = SHPdataStore.getFeatureSource();
    ContentFeatureCollection featureCollection = featureSource.getFeatures();
    SimpleFeatureType schema = featureSource.getSchema();
    CoordinateReferenceSystem sourceCRS = schema.getCoordinateReferenceSystem();
    System.out.println(sourceCRS);
    // the coordinates system I want to reproject the data to
    CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
    // Define a MathTransform object
    MathTransform transform = CRS.findMathTransform(sourceCRS, targetCRS);
    try (SimpleFeatureIterator iterator = featureCollection.features()) {
      while (iterator.hasNext()) {
        SimpleFeature type = iterator.next();
        Geometry defaultGeometry = (Geometry) type.getDefaultGeometry();
        System.out.println(defaultGeometry);
        System.out.println(JTS.transform(defaultGeometry, transform));
      }
    }

  }

}

Output

ShapefileDataStore [file=https://github.com/basilmohammed/CRS-Extraction-From-ShapeFile/raw/master/ShapeFiles/3857/3857-Point/3857Point.shp, charset=ISO-8859-1, timeZone=sun.util.calendar.ZoneInfo[id="Europe/London",offset=0,dstSavings=3600000,useDaylight=true,transitions=242,lastRule=java.util.SimpleTimeZone[id=Europe/London,offset=0,dstSavings=3600000,useDaylight=true,startYear=0,startMode=2,startMonth=2,startDay=-1,startDayOfWeek=1,startTime=3600000,startTimeMode=2,endMode=2,endMonth=9,endDay=-1,endDayOfWeek=1,endTime=3600000,endTimeMode=2]], memoryMapped=false, bufferCachingEnabled=true, indexed=true, fidIndexed=true]
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere", 
  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["Popular Visualisation Pseudo Mercator"], 
  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["x", EAST], 
  AXIS["y", NORTH]]
POINT (77.39083006594812 23.27470458279848)
POINT (0.0002090802286007 0.0006952136549894)
...

Thus the first point (POINT (77.39083006594812 23.27470458279848)) is transformed to POINT (0.0002090802286007 0.0006952136549894) which is probably not what you want.

If I switch to .../ShapeFiles/4326/4326-Point/point.shp and transform to EPSG:3857 I get: POINT (77.39083006594812 23.27470458279848) -> POINT (8615107.79501011 2665273.4881518427) which looks plausible.