-1

I have been trying to change the projection of Shapefile from one coordinate reference system to other. The shapefile I have used has EPSG:4326 as its reference system and I need to change it to EPSG:32056. I am using Geotools API-20.0 for the same.

I have already tried various methods available in the geotools like using ReprojectingFeatureCollection, use of JTS, use of Query API to convert the shapefile directly to the other coordinate reference system

   import java.awt.event.ActionEvent;
   import java.io.File;
   import java.io.InputStream;
   import java.io.Serializable;
   import java.net.URI;
   import java.net.URL;

    import java.util.HashMap;
    import java.util.Map;
    import java.util.Set;

    import javax.swing.Action;
    import javax.swing.JButton;
    import javax.swing.JOptionPane;
    import javax.swing.JToolBar;
    import javax.swing.SwingWorker;

    import org.geotools.data.DataStore;

import org.geotools.data.DataStoreFactorySpi;
import org.geotools.data.DefaultTransaction;
import org.geotools.data.FeatureWriter;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.data.Query;
import org.geotools.data.Transaction;
import org.geotools.data.shapefile.ShapefileDataStoreFactory;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.data.simple.SimpleFeatureStore;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geometry.jts.JTS;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.map.FeatureLayer;
import org.geotools.map.Layer;
import org.geotools.map.MapContent;
import org.geotools.referencing.CRS;
import org.geotools.referencing.ReferencingFactoryFinder;
import org.geotools.referencing.factory.gridshift.GridShiftLocator;
import org.geotools.styling.SLD;
import org.geotools.styling.Style;
import org.geotools.swing.JMapFrame;
import org.geotools.swing.JProgressWindow;
import org.geotools.swing.action.SafeAction;
import org.geotools.swing.data.JFileDataStoreChooser;
import org.locationtech.jts.geom.Envelope;
import org.opengis.feature.Feature;
import org.opengis.feature.FeatureVisitor;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.feature.type.FeatureType;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.util.ProgressListener;

import com.vividsolutions.jts.geom.Geometry;

public class CRSLab {

    private File sourceFile;
    private SimpleFeatureSource featureSource;
    private MapContent map;

    public static void main(String[] args) throws Exception {
        CRSLab lab = new CRSLab();
        lab.displayShapefile();
    }

    // docs end main
    /**
     * This method:
     * <ol type="1">
     * <li>Prompts the user for a shapefile to display
     * <li>Creates a JMapFrame with custom toolbar buttons
     * <li>Displays the shapefile
     * </ol>
     */
    // docs start display
    private void displayShapefile() throws Exception {
        sourceFile = JFileDataStoreChooser.showOpenFile("shp", null);
        if (sourceFile == null) {
            return;
        }
        FileDataStore store = FileDataStoreFinder.getDataStore(sourceFile);
        featureSource = store.getFeatureSource();

        // Create a map context and add our shapefile to it
        map = new MapContent();
        Style style = SLD.createSimpleStyle(featureSource.getSchema());
        Layer layer = new FeatureLayer(featureSource, style);
        map.layers().add(layer);

        // Create a JMapFrame with custom toolbar buttons
        JMapFrame mapFrame = new JMapFrame(map);
        mapFrame.enableToolBar(true);
        mapFrame.enableStatusBar(true);

        JToolBar toolbar = mapFrame.getToolBar();
        toolbar.addSeparator();
        toolbar.add(new JButton(new ValidateGeometryAction()));
        toolbar.add(new JButton(new ExportShapefileAction()));

        // Display the map frame. When it is closed the application will exit
        mapFrame.setSize(800, 600);
        mapFrame.setVisible(true);
    }
    // docs end display


    // docs start export
    private void exportToShapefile() throws Exception {
        SimpleFeatureType schema = featureSource.getSchema();
        JFileDataStoreChooser chooser = new JFileDataStoreChooser("shp");
        chooser.setDialogTitle("Save reprojected shapefile");
        chooser.setSaveFile(sourceFile);
        int returnVal = chooser.showSaveDialog(null);
        if (returnVal != JFileDataStoreChooser.APPROVE_OPTION) {
            return;
        }
        File file = chooser.getSelectedFile();
        if (file.equals(sourceFile)) {
            JOptionPane.showMessageDialog(null, "Cannot replace " + file);
            return;
        }

        // set up the math transform used to process the data
        CoordinateReferenceSystem dataCRS = schema.getCoordinateReferenceSystem();
        CoordinateReferenceSystem worldCRS = CRS.decode("EPSG:32056", true);// map.getCoordinateReferenceSystem();
        boolean lenient = true; // allow for some error due to different datums
        MathTransform transform = CRS.findMathTransform(dataCRS, worldCRS, lenient);

        // grab all features
        SimpleFeatureCollection featureCollection = featureSource.getFeatures();

        // And create a new Shapefile with a slight modified schema
        DataStoreFactorySpi factory = new ShapefileDataStoreFactory();
        Map<String, Serializable> create = new HashMap<String, Serializable>();
        create.put("url", file.toURI().toURL());
        create.put("create spatial index", Boolean.TRUE);
        DataStore dataStore = factory.createNewDataStore(create);
        SimpleFeatureType featureType = SimpleFeatureTypeBuilder.retype(schema, worldCRS);
        dataStore.createSchema(featureType);
        String createdName = dataStore.getTypeNames()[0];
        // carefully open an iterator and writer to process the results
        Transaction transaction = new DefaultTransaction("Reproject");
        FeatureWriter<SimpleFeatureType, SimpleFeature> writer = dataStore.getFeatureWriterAppend(createdName,
                transaction);
        SimpleFeatureIterator iterator = featureCollection.features();
        try {
            int counter = 0;
            while (iterator.hasNext()) {
                // copy the contents of each feature and transform the geometry
                SimpleFeature feature = iterator.next();
                SimpleFeature copy = writer.next();

                org.locationtech.jts.geom.Geometry geometry = (org.locationtech.jts.geom.Geometry) feature
                        .getDefaultGeometry();
                org.locationtech.jts.geom.Geometry geometry2 = JTS.transform(geometry, transform);
                System.out.println(geometry.isSimple() && geometry2.isSimple());
                // if (geometry2.isValid()) {
                copy.setAttributes(feature.getAttributes());
                counter++;
                copy.setDefaultGeometry(geometry2);
                writer.write();
                // }
            }
            transaction.commit();
            System.out.println("valid geometries : " + counter);
            JOptionPane.showMessageDialog(null, "Export to shapefile complete");
        } catch (Exception problem) {
            problem.printStackTrace();
            transaction.rollback();
            JOptionPane.showMessageDialog(null, "Export to shapefile failed");
        } finally {
            writer.close();
            iterator.close();
            transaction.close();
        }
    }
    // docs end export

        // docs start validate
    private int validateFeatureGeometry(ProgressListener progress) throws Exception {
        final SimpleFeatureCollection featureCollection = featureSource.getFeatures();

        // Rather than use an iterator, create a FeatureVisitor to check each
        // fature
        class ValidationVisitor implements FeatureVisitor {
            public int numInvalidGeometries = 0;

            public void visit(Feature f) {
                SimpleFeature feature = (SimpleFeature) f;
                Geometry geom = (Geometry) feature.getDefaultGeometry();
                if (geom != null && !geom.isValid()) {
                    numInvalidGeometries++;
                    System.out.println("Invalid Geoemtry: " + feature.getID());
                }
            }
        }

        ValidationVisitor visitor = new ValidationVisitor();

        // Pass visitor and the progress bar to feature collection
        featureCollection.accepts(visitor, progress);
        return visitor.numInvalidGeometries;
    }
    // docs end validate

        // docs start export action
    class ExportShapefileAction extends SafeAction {
        ExportShapefileAction() {
            super("Export...");
            putValue(Action.SHORT_DESCRIPTION, "Export using current crs");
        }

        public void action(ActionEvent e) throws Throwable {
            exportToShapefile();
        }
    }
    // docs end export action

    /**
     * This class performs the task of checking that the Geometry of each
     * feature is topologically valid and reports on the results. It also
     * supplies the name and tool tip.
     */
    // docs start validate action
    class ValidateGeometryAction extends SafeAction {
        ValidateGeometryAction() {
            super("Validate geometry");
            putValue(Action.SHORT_DESCRIPTION, "Check each geometry");
        }

        public void action(ActionEvent e) throws Throwable {
            int numInvalid = validateFeatureGeometry(null);
            String msg;
            if (numInvalid == 0) {
                msg = "All feature geometries are valid";
            } else {
                msg = "Invalid geometries: " + numInvalid;
            }
            JOptionPane.showMessageDialog(null, msg, "Geometry results", JOptionPane.INFORMATION_MESSAGE);
        }
    }
        // docs end validate action

    }

The output obtained after doing projection using Geotools are a lot different than what I used to get from ArcMap of esri. Is there any other transformation that I should perform.

GeoTools_Output

ArcMapOutput

Ian Turton
  • 10,018
  • 1
  • 28
  • 47
  • can you show a "correct" image too, or provide a link to the input shapefile? – Ian Turton Apr 20 '19 at 07:45
  • @IanTurton: I have uploaded the expected output as ArcMapOutput and the shapefile i used can be downloaded from https://hub.arcgis.com/datasets/a21fdb46d23e4ef896f31475217cbb08_1 – naveen nagar Apr 23 '19 at 06:13

1 Answers1

1

When I try this (with v22.x) all I get is an error as too many points are outside the valid projection error. This is because you are taking a map of the world and reprojecting it to a CRS designed for Wyoming.

It seems that ESRI are being "helpful" and clipping your output to the area of validity (assuming you meant something other than EPSG:32056). GeoTools assumes that you know what you are doing and doesn't do that, which is why you have all the countries of the world shown in that map.

Here is the output for just the USA, which suggests that the ESRI image you show is a different projection again (look at the 49th parallel).

enter image description here

Ian Turton
  • 10,018
  • 1
  • 28
  • 47
  • Hi Ian, I am new to the Geotools, Could you please let me know if there is any way where we can get the valid extent of the projection system. As I will be requiring to project to other projections as well or shall I use isValid function at org.locationtech.jts.geom.Geometry in order to only show the correct data which lies in the valid extent. – naveen nagar Apr 24 '19 at 15:32