Skip to content

Instantly share code, notes, and snippets.

@ianturton
Created November 19, 2020 11:43
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save ianturton/aec5ddc431de156394f0a131f1a61ad0 to your computer and use it in GitHub Desktop.
Save ianturton/aec5ddc431de156394f0a131f1a61ad0 to your computer and use it in GitHub Desktop.
Simple GeoTools code to generate the Null Island archipelago
package com.ianturton.cookbook.projections;
import java.io.File;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Collection;
import java.util.HashSet;
import java.util.Set;
import org.apache.commons.cli.CommandLine;
import org.apache.commons.cli.CommandLineParser;
import org.apache.commons.cli.DefaultParser;
import org.apache.commons.cli.Options;
import org.apache.commons.cli.ParseException;
import org.geotools.data.DataUtilities;
import org.geotools.data.shapefile.ShapefileDataStore;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureStore;
import org.geotools.feature.FeatureCollection;
import org.geotools.feature.simple.SimpleFeatureBuilder;
import org.geotools.feature.simple.SimpleFeatureTypeBuilder;
import org.geotools.geometry.jts.JTS;
import org.geotools.metadata.iso.IdentifierImpl;
import org.geotools.referencing.CRS;
import org.geotools.referencing.ReferencingFactoryFinder;
import org.geotools.referencing.crs.DefaultGeographicCRS;
import org.geotools.util.URLs;
import org.geotools.util.factory.Hints;
import org.locationtech.jts.geom.Coordinate;
import org.locationtech.jts.geom.GeometryFactory;
import org.locationtech.jts.geom.Point;
import org.opengis.feature.simple.SimpleFeature;
import org.opengis.feature.simple.SimpleFeatureType;
import org.opengis.metadata.Identifier;
import org.opengis.referencing.FactoryException;
import org.opengis.referencing.NoSuchAuthorityCodeException;
import org.opengis.referencing.crs.CRSAuthorityFactory;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.crs.ProjectedCRS;
import org.opengis.referencing.operation.MathTransform;
/**
* Find all of the Null Islands in the World cf https://t.co/OZ9s20LQAx
*
* @author ian
*
*/
public class NullIslands {
private static FeatureCollection featureCollection;
public static void main(String[] args) throws NoSuchAuthorityCodeException, FactoryException, IOException {
NullIslands guesser = new NullIslands();
CommandLineParser parser = new DefaultParser();
Options options = new Options();
options.addOption("e", true, "Limit search to just this EPSG code");
CommandLine commandLine;
try {
commandLine = parser.parse(options, args);
if (commandLine.getOptions().length == 0) {
guesser.setAll(true);
}
if (commandLine.hasOption('e')) {
guesser.setAll(false);
guesser.setCRS(CRS.decode(commandLine.getOptionValue('e')));
}
SimpleFeatureCollection results = guesser.search();
File output = new File("nullislands.shp");
ShapefileDataStore ds = new ShapefileDataStore(URLs.fileToUrl(output));
ds.createSchema(guesser.getType());
SimpleFeatureStore store = (SimpleFeatureStore) ds.getFeatureSource();
store.addFeatures(results);
} catch (ParseException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
private SimpleFeatureType type;
private SimpleFeatureType getType() {
if (this.type == null) {
SimpleFeatureTypeBuilder ftBuilder = new SimpleFeatureTypeBuilder();
ftBuilder.setName("NullPoints");
ftBuilder.add("code", String.class);
ftBuilder.add("the_geom", Point.class);
ftBuilder.setCRS(DefaultGeographicCRS.WGS84);
ftBuilder.setDefaultGeometry("the_geom");
type = ftBuilder.buildFeatureType();
}
return type;
}
private CoordinateReferenceSystem crs;
private void setCRS(CoordinateReferenceSystem decode) {
this.crs = decode;
}
private boolean all = true;
private void setAll(boolean b) {
this.all = b;
}
public SimpleFeatureCollection search() {
//
//
ArrayList<SimpleFeature> ret = new ArrayList<>();
Set<ProjectedCRS> possibles = new HashSet<>();
GeometryFactory geometryFactory = new GeometryFactory();
Point nullPoint = geometryFactory.createPoint(new Coordinate(0, 0));
SimpleFeatureBuilder builder = new SimpleFeatureBuilder(getType());
if (isAll()) {
Set<String> auths = ReferencingFactoryFinder.getAuthorityNames();
for (String auth : auths) {
possibles.addAll(fetchPossibles(auth));
}
} else {
possibles.add((ProjectedCRS) this.crs);
}
for (ProjectedCRS proj : possibles) {
try {
MathTransform testTransform = CRS.findMathTransform(proj, DefaultGeographicCRS.WGS84);
Point p = (Point) JTS.transform(nullPoint, testTransform);
builder.set("code", CRS.lookupIdentifier(proj, false));
builder.set("the_geom", p);
ret.add(builder.buildFeature(null));
} catch (Exception e) {
// lots of things can go wrong with random tests so ignore them
}
}
return DataUtilities.collection(ret);
}
/**
* @return
*/
private Set<ProjectedCRS> fetchPossibles(String authority) {
Set<ProjectedCRS> possibles = new HashSet<>();
Hints hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.FALSE);
Set<CRSAuthorityFactory> facs = ReferencingFactoryFinder.getCRSAuthorityFactories(hints);
Set<String> codes = new HashSet<>();
CRSAuthorityFactory factory = null;
int bestCodes = 0;
Set<String> c;
IdentifierImpl identifier = new IdentifierImpl(authority);
for (CRSAuthorityFactory fac : facs) {
Collection<? extends Identifier> identifiers = fac.getAuthority().getIdentifiers();
if (identifiers.contains(identifier)) {
try {
c = fac.getAuthorityCodes(org.opengis.referencing.crs.ProjectedCRS.class);
if (c.size() > bestCodes) {
bestCodes = c.size();
factory = fac;
codes = c;
}
} catch (FactoryException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
}
for (String code : codes) {
try {
ProjectedCRS crs = factory.createProjectedCRS(code);
possibles.add(crs);
} catch (FactoryException e) {
// TODO Auto-generated catch block
// e.printStackTrace();
}
}
return possibles;
}
public boolean isAll() {
return all;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment