1

我正在尝试编写一个 Java 程序,它将大量 GPS 坐标捕捉到线 shapefile(道路网络),并且不仅返回新坐标,还返回捕捉到的线段的唯一标识符。该标识符是 FID、其他语言中使用的“索引”(即 1 是第一个特征等)或属性表中的任何列都没有关系。

我已经在 R 中使用maptools::snapPointsToLines函数完成了这项工作,但是考虑到我需要处理的数据量,这是不可扩展的,所以我正在寻找 Java 来更快地处理数据以便在 R 中进行分析。

我的代码(如下)目前与用于捕捉的 geotools 教程非常相似,只是我在 GPS 点的(1900 万行)CSV 中读取而不是生成它们,并且我编写了结果的 CSV。它捕捉得很好,而且比我得到的要快得多,但我不知道如何识别捕捉到的线。可用的文档似乎涵盖了对功能集的查询和过滤,我无法使其特别适用于该代码创建的索引行对象,并且我的代码中的现有函数toString()返回了一些我无法理解的内容,例如com.vividsolutions.jts.linearreff.LocationIndexedLine@74cec793.

基本上,我只想让 lineID 字段生成任何其他 GIS 软件或语言都可以匹配特定路段的内容。

package org.geotools.tutorial.quickstart;

import java.io.*;
import java.util.List;
import java.util.Arrays;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Envelope;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.MultiLineString;
import com.vividsolutions.jts.index.SpatialIndex;
import com.vividsolutions.jts.index.strtree.STRtree;
import com.vividsolutions.jts.linearref.LinearLocation;
import com.vividsolutions.jts.linearref.LocationIndexedLine;

import org.geotools.data.FeatureSource;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.feature.FeatureCollection;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.swing.data.JFileDataStoreChooser;
import org.geotools.util.NullProgressListener;
import org.opengis.feature.Feature;
import org.opengis.feature.FeatureVisitor;
import org.opengis.feature.simple.SimpleFeature;
import com.opencsv.*;

public class SnapToLine {

    public static void main(String[] args) throws Exception {

        /*
         * Open a shapefile. You should choose one with line features
         * (LineString or MultiLineString geometry)
         * 
         */
        File file = JFileDataStoreChooser.showOpenFile("shp", null);
        if (file == null) {
            return;
        }

        FileDataStore store = FileDataStoreFinder.getDataStore(file);
        FeatureSource source = store.getFeatureSource();

        // Check that we have line features
        Class<?> geomBinding = source.getSchema().getGeometryDescriptor().getType().getBinding();
        boolean isLine = geomBinding != null 
                && (LineString.class.isAssignableFrom(geomBinding) ||
                    MultiLineString.class.isAssignableFrom(geomBinding));

        if (!isLine) {
            System.out.println("This example needs a shapefile with line features");
            return;
        }
         final SpatialIndex index = new STRtree();
        FeatureCollection features = source.getFeatures();
        //FeatureCollection featurecollection = source.getFeatures(Query.FIDS);
        System.out.println("Slurping in features ...");
        features.accepts(new FeatureVisitor() {

            @Override
            public void visit(Feature feature) {
                SimpleFeature simpleFeature = (SimpleFeature) feature;
                Geometry geom = (MultiLineString) simpleFeature.getDefaultGeometry();
                // Just in case: check for  null or empty geometry
                if (geom != null) {
                    Envelope env = geom.getEnvelopeInternal();
                    if (!env.isNull()) {
                        index.insert(env, new LocationIndexedLine(geom));
                    }
                }
            }
        }, new NullProgressListener());
 /*

 /*
         * We defined the maximum distance that a line can be from a point
         * to be a candidate for snapping 
         */

        ReferencedEnvelope bounds = features.getBounds();
        final double MAX_SEARCH_DISTANCE = bounds.getSpan(0) / 1000.0;



        int pointsProcessed = 0;
        int pointsSnapped = 0;
        long elapsedTime = 0;
        long startTime = System.currentTimeMillis();
        double longiOut;
        double latiOut;
        int moved;
        String lineID   = "NA";

        //Open up the CSVReader. Reading in line by line to avoid memory failure.

        CSVReader csvReader = new CSVReader(new FileReader(new File("fakedata.csv")));
        String[] rowIn;



        //open up the CSVwriter
        String outcsv = "fakedataOUT.csv";
        CSVWriter writer = new CSVWriter(new FileWriter(outcsv));



        while ((rowIn = csvReader.readNext()) != null) {

            // Get point and create search envelope
            pointsProcessed++;
            double longi = Double.parseDouble(rowIn[0]);
            double lati  = Double.parseDouble(rowIn[1]);
            Coordinate pt = new Coordinate(longi, lati);
            Envelope search = new Envelope(pt);
            search.expandBy(MAX_SEARCH_DISTANCE);

            /*
             * Query the spatial index for objects within the search envelope.
             * Note that this just compares the point envelope to the line envelopes
             * so it is possible that the point is actually more distant than
             * MAX_SEARCH_DISTANCE from a line.
             */
            List<LocationIndexedLine> lines = index.query(search);

            // Initialize the minimum distance found to our maximum acceptable
            // distance plus a little bit
            double minDist = MAX_SEARCH_DISTANCE + 1.0e-6;
            Coordinate minDistPoint = null;

            for (LocationIndexedLine line : lines) {
                LinearLocation here = line.project(pt);
                Coordinate point = line.extractPoint(here);
                double dist = point.distance(pt);
                if (dist < minDist) {
                    minDist = dist;
                    minDistPoint = point;
                    lineID = line.toString();
                }
            }


            if (minDistPoint == null) {
                // No line close enough to snap the point to
                System.out.println(pt + "- X");
                longiOut = longi;
                latiOut  = lati;
                moved    = 0;
                lineID   = "NA";
            } else {
                System.out.printf("%s - snapped by moving %.4f\n", 
                        pt.toString(), minDist);
                longiOut = minDistPoint.x;
                latiOut  = minDistPoint.y;
                moved    = 1;        
                pointsSnapped++;
            }
    //write a new row

    String [] rowOut = {Double.toString(longiOut), Double.toString(latiOut), Integer.toString(moved), lineID}; 
    writer.writeNext(rowOut);
        }

        System.out.printf("Processed %d points (%.2f points per second). \n"
                + "Snapped %d points.\n\n",
                pointsProcessed,
                1000.0 * pointsProcessed / elapsedTime,
                pointsSnapped);
        writer.close();
    }
}

我不仅是 Java 新手,而且只是在 R 等特定领域的语言中自学过;我不是一个编码器,而是一个使用代码的人,所以如果解决方案看起来很明显,我可能缺乏基本理论!

ps我知道那里有更好的地图匹配解决方案(graphhopper等),我只是想开始eas!

谢谢!

4

2 回答 2

0

我会尽量避免深入 JTS 兔子洞并坚持使用 GeoTools(当然我是 GeoTools 开发人员,所以我会这么说)。

首先,我会使用 anSpatialIndexFeatureCollection来保存我的行(假设它们适合内存,否则 PostGIS 表是要走的路)。这使我不必建立自己的索引。

然后我会用 aCSVDataStore来保存从 GPS 流中解析我自己的点(因为我很懒,那里也有很多问题)。

这意味着大部分工作归结为这个循环,DWITHIN找到指定距离内的所有特征:

try (SimpleFeatureIterator itr = pointFeatures.getFeatures().features()) { 
  while (itr.hasNext()) {
    SimpleFeature f = itr.next();
    Geometry snapee = (Geometry) f.getDefaultGeometry();
    Filter filter = ECQL.toFilter("DWITH(\"the_geom\",'" + writer.write(snapee) + "'," + MAX_SEARCH_DISTANCE + ")");
    SimpleFeatureCollection possibles = indexed.subCollection(filter);
    double minDist = Double.POSITIVE_INFINITY;
    SimpleFeature bestFit = null;
    Coordinate bestPoint = null;
    try (SimpleFeatureIterator pItr = possibles.features()) {
      while (pItr.hasNext()) {
        SimpleFeature p = pItr.next();
        Geometry line = (Geometry) p.getDefaultGeometry();

        double dist = snapee.distance(line);
        if (dist < minDist) {
          minDist = dist;
          bestPoint = DistanceOp.nearestPoints(snapee, line)[1];
          bestFit  = p;
        }
      }
    }

在该循环结束时,您应该从线条(包括其 id 和名称等)、最近点 (bestPoint) 和移动距离 (minDist) 中知道最近的特征 (bestFit)。

同样,我可能也会使用 aCSVDatastore将功能写回。

如果您有数百万个点,我可能会考虑使用 aFilterFactory直接创建过滤器,而不是使用ECQL解析器。

于 2017-01-11T18:01:25.217 回答
0

根据 iant 在接受的答案中提供的代码,我已将代码更改如下,我将其发布给其他任何在谷歌上搜索此类问题的人。

几件事:确保ECQL解析器中有一个单元(我还没有尝试过这个FilterFactory建议)。以下CSVDataStore链接中介绍了这一点。请注意,默认情况下,经度和纬度分别硬编码为lonlathttp://docs.geotools.org/latest/userguide/tutorial/datastore/intro.html

package org.geotools.tutorial.quickstart;

import java.io.*;
import java.util.List;
import java.util.Arrays;

import com.vividsolutions.jts.geom.Coordinate;
import com.vividsolutions.jts.geom.Envelope;
import com.vividsolutions.jts.geom.Geometry;
import com.vividsolutions.jts.geom.LineString;
import com.vividsolutions.jts.geom.MultiLineString;
import com.vividsolutions.jts.index.SpatialIndex;
import com.vividsolutions.jts.index.strtree.STRtree;
import com.vividsolutions.jts.linearref.LinearLocation;
import com.vividsolutions.jts.linearref.LocationIndexedLine;
import com.vividsolutions.jts.operation.distance.DistanceOp;
import com.vividsolutions.jts.io.WKTWriter;

import org.geotools.data.FeatureSource;
import org.geotools.data.simple.SimpleFeatureSource;
import org.geotools.data.simple.SimpleFeatureCollection;
import org.geotools.data.simple.SimpleFeatureIterator;
import org.geotools.data.FileDataStore;
import org.geotools.data.FileDataStoreFinder;
import org.geotools.feature.FeatureCollection;
import org.geotools.geometry.jts.ReferencedEnvelope;
import org.geotools.swing.data.JFileDataStoreChooser;
import org.geotools.util.NullProgressListener;
import org.geotools.data.collection.SpatialIndexFeatureCollection;
import org.geotools.data.collection.SpatialIndexFeatureSource;
import org.geotools.filter.text.ecql.ECQL;

import org.opengis.filter.Filter;
import org.opengis.feature.Feature;
import org.opengis.feature.FeatureVisitor;
import org.opengis.feature.simple.SimpleFeature;
import com.opencsv.*;
import com.csvreader.CsvReader;




        public class SnapToLine {

    public static void main(String[] args) throws Exception {

        //input and output files and other parameters
        String inputpoints = "/home/bitre/fakedata.csv";
        String outcsv = "fakedataOUT.csv";

        final double MAX_SEARCH_DISTANCE = 0.5;


        /*
         * Open a shapefile. You should choose one with line features
         * (LineString or MultiLineString geometry)
         * 
         */
        File file = JFileDataStoreChooser.showOpenFile("shp", null);
        if (file == null) {
            return;
        }

        FileDataStore store = FileDataStoreFinder.getDataStore(file);
        SimpleFeatureSource source = store.getFeatureSource();

        // Check that we have line features
        Class<?> geomBinding = source.getSchema().getGeometryDescriptor().getType().getBinding();
        boolean isLine = geomBinding != null 
                && (LineString.class.isAssignableFrom(geomBinding) ||
                    MultiLineString.class.isAssignableFrom(geomBinding));

        if (!isLine) {
            System.out.println("This example needs a shapefile with line features");
            return;
        }



        SimpleFeatureCollection features = source.getFeatures();

        SpatialIndexFeatureCollection indexed = new SpatialIndexFeatureCollection(features);



 /*

 /*
         * We defined the maximum distance that a line can be from a point
         * to be a candidate for snapping 
         */

        ReferencedEnvelope bounds = features.getBounds();




        //open up the CSVwriter

        CSVWriter csvWriter = new CSVWriter(new FileWriter(outcsv));



        //CSVDataStore features for the points


        CSVDataStore pointFeaturesCSV = new CSVDataStore(new File(inputpoints));
        String typeName = pointFeaturesCSV.getTypeNames()[0];
        SimpleFeatureSource pointFeatures = pointFeaturesCSV.getFeatureSource(typeName); 


        double longiOut;
        double latiOut;
        int progress = 0;
        int remn;
        String[] rowOut = new String[4];

                try (SimpleFeatureIterator itr = pointFeatures.getFeatures().features()) { 
           while (itr.hasNext()) {
             SimpleFeature f = itr.next();
             Geometry snapee = (Geometry) f.getDefaultGeometry();
             WKTWriter writer = new WKTWriter();
             Filter filter = ECQL.toFilter("DWITHIN(\"the_geom\",'" + writer.write(snapee) + "'," + MAX_SEARCH_DISTANCE + "," + "kilometers" + ")");
             SimpleFeatureCollection possibles = indexed.subCollection(filter);
             double minDist = Double.POSITIVE_INFINITY;
             SimpleFeature bestFit = null;
             Coordinate bestPoint = null;
             try (SimpleFeatureIterator pItr = possibles.features()) {
               while (pItr.hasNext()) {
                 SimpleFeature p = pItr.next();
                 Geometry line = (Geometry) p.getDefaultGeometry();

                 double dist = snapee.distance(line);
                 if (dist < minDist) {
                   minDist = dist;
                   bestPoint = DistanceOp.nearestPoints(snapee, line)[1]; // google DistanceOp
                   bestFit  = p;
                 }
                 longiOut = bestPoint.x;
                 latiOut = bestPoint.y;
                 rowOut[0] = bestFit.getID();
                 rowOut[1] = Double.toString(minDist);
                 rowOut[2] = Double.toString(longiOut);
                 rowOut[3] = Double.toString(latiOut);

                 //rowOut = {bestFit.getID(), Double.toString(minDist), Double.toString(longiOut), Double.toString(latiOut)};

               }
                 csvWriter.writeNext(rowOut);
                progress ++;
                remn = progress % 1000000;
                if(remn == 0){
                    System.out.println("Just snapped line" + progress);
                }

             }

           }


       }
    }
}
于 2017-01-16T03:46:52.797 回答