|
1 | 1 | package org.geoscript.example
|
2 | 2 |
|
3 |
| -import org.geoscript._ |
4 |
| - |
5 |
| -object Intersections { |
6 |
| - def process(src: layer.Layer, dest: layer.Layer, joinField: String) { |
7 |
| - println("Processing %s".format(src.schema.name)) |
8 |
| - |
9 |
| - for (feat <- src.features) { |
10 |
| - val intersections = |
11 |
| - src.filter(filter.Filter.intersects(feat.geometry)) |
12 |
| - dest ++= |
13 |
| - intersections.filter(_.id > feat.id).map { corner => |
14 |
| - feature.Feature( |
15 |
| - "geom" -> (feat.geometry intersection corner.geometry), |
16 |
| - (joinField + "Left") -> feat.get[Any](joinField), |
17 |
| - (joinField + "Right") -> corner.get[Any](joinField) |
18 |
| - ) |
19 |
| - } |
20 |
| - } |
21 |
| - |
22 |
| - println("Found %d intersections".format(dest.count)) |
23 |
| - } |
| 3 | +import org.geoscript.feature._ |
| 4 | +import org.geoscript.filter._ |
| 5 | +import org.geoscript.filter.builder._ |
| 6 | +import org.geoscript.geometry._ |
| 7 | +import org.geoscript.layer._ |
| 8 | +import org.geoscript.workspace._ |
| 9 | + |
| 10 | +/** |
| 11 | + * An example app for creating a shapefile containing all intersections between |
| 12 | + * two input Shapefiles |
| 13 | + */ |
| 14 | +object Intersections extends App { |
| 15 | + // by extending the scala.App trait from the Scala standard library, we avoid |
| 16 | + // writing the def main(args: Array[String]) that is otherwise required for |
| 17 | + // an object to be executable. |
24 | 18 |
|
25 |
| - def rewrite(schema: feature.Schema, fieldName: String): feature.Schema = |
26 |
| - feature.Schema( |
27 |
| - schema.name + "_intersections", |
28 |
| - feature.Field( |
29 |
| - "geom", |
30 |
| - classOf[com.vividsolutions.jts.geom.Geometry], |
31 |
| - schema.geometry.projection |
32 |
| - ), |
33 |
| - feature.Field(fieldName + "Left", classOf[String]), |
34 |
| - feature.Field(fieldName + "Right", classOf[String]) |
35 |
| - ) |
36 |
| - |
37 |
| - def main(args: Array[String]) = { |
38 |
| - if (args.length == 0) { |
39 |
| - println("You need to provide the path to a shapefile as an argument to this example.") |
40 |
| - } else { |
41 |
| - val src = layer.Shapefile(args(0)) |
42 |
| - val joinField = |
43 |
| - src.schema.fields.find { _.gtBinding == classOf[String] } match { |
44 |
| - case Some(f) => f.name |
45 |
| - case None => "id" |
46 |
| - } |
47 |
| - val dest = src.workspace.create(rewrite(src.schema, joinField)) |
48 |
| - process(src, dest, joinField) |
49 |
| - } |
| 19 | + if (args.size < 3) { |
| 20 | + // We need three arguments: first shapefile to scan, second shapefile to |
| 21 | + // scan, shapefile to store results. |
| 22 | + println( |
| 23 | +""" |Usage: Intersections <first shapefile> <second shapefile> <output shapefile> |
| 24 | + |Computes all intersections between the two input shapefiles and stores the |
| 25 | + |results in the output shapefile. The output will also have two fields |
| 26 | + |named left_id and right_id containing the ids of the features that |
| 27 | + |intersected. (This is just an example - NOTE that shapefile features do not |
| 28 | + |have stable identifiers.)""".stripMargin) |
| 29 | + System.exit(0) |
50 | 30 | }
|
| 31 | + |
| 32 | + // for convenience, we create a little function for connecting to shapefiles. |
| 33 | + val connect = (path: String) => |
| 34 | + new org.geotools.data.shapefile.ShapefileDataStore( |
| 35 | + new java.io.File(path).toURI.toURL): Workspace |
| 36 | + |
| 37 | + // Here we use a pattern match to concisely extract the arguments into |
| 38 | + // individual variables. |
| 39 | + val Array(leftPath, rightPath, outputPath) = (args take 3) |
| 40 | + val leftLayer = connect(leftPath).layers.head |
| 41 | + val rightLayer = connect(rightPath).layers.head |
| 42 | + |
| 43 | + // There are a few different ways to compute the intersections. |
| 44 | + // The simplest is to use a Scala for-comprehension. |
| 45 | + // val intersections = |
| 46 | + // for { |
| 47 | + // l <- leftLayer.features |
| 48 | + // r <- rightLayer.features |
| 49 | + // if l.geometry intersects r.geometry |
| 50 | + // } yield (l.geometry intersection r.geometry, l.id, r.id) |
| 51 | + |
| 52 | + // This produces correct results, but there are some performance problems. |
| 53 | + // * It fetches all features from the 'right' layer on each step of iterating |
| 54 | + // through the 'left' layer. This might mean a lot of disk access! |
| 55 | + // * The results are stored in memory. Since we're just going to write the |
| 56 | + // features to a new shapefile it would be nice to avoid that. It would |
| 57 | + // save some memory, and also might complete faster if we can start writing |
| 58 | + // the results before we finish finding all the intersections. |
| 59 | + // |
| 60 | + // We can avoid repetitive requests to the underlying store by copying all the |
| 61 | + // features into an in-memory collection before scanning. |
| 62 | + |
| 63 | + // val leftFeatures = leftLayer.features.to[Vector] |
| 64 | + // val rightFeatures = rightLayer.features.to[Vector] |
| 65 | + |
| 66 | + // val intersections2 = |
| 67 | + // for { |
| 68 | + // l <- leftFeatures |
| 69 | + // r <- rightFeatures |
| 70 | + // if l.geometry intersects r.geometry |
| 71 | + // } yield (l.geometry intersection r.geometry, l.id, r.id) |
| 72 | + |
| 73 | + // This trades off memory in order to speed up the processing, so it's |
| 74 | + // still only going to work for small datasets. Instead of performing the |
| 75 | + // filtering in Scala code, we can use the GeoTools Query system to have it |
| 76 | + // performed by the datastore itself. Depending on the datastore filters will |
| 77 | + // be more or less completely executed by the underlying engine. For example, |
| 78 | + // filters executed against a Postgis database can be largely converted to |
| 79 | + // SQL. For Shapefiles most filter operations are executed in-process, but |
| 80 | + // they are at least able to take advantage of a spatial index. |
| 81 | + |
| 82 | + val intersections3 = |
| 83 | + for { |
| 84 | + l <- leftLayer.features |
| 85 | + r <- rightLayer.filter( |
| 86 | + Literal(l.geometry) intersects Property(rightLayer.schema.geometryField.name)) |
| 87 | + } yield (l.geometry intersection r.geometry, l.id, r.id) |
| 88 | + |
| 89 | + intersections3.foreach(println) |
| 90 | + |
| 91 | + // require(intersections.toSeq == intersections.toSeq.distinct) |
| 92 | + |
| 93 | + // def process(src: layer.Layer, dest: layer.Layer, joinField: String) { |
| 94 | + // println("Processing %s".format(src.schema.name)) |
| 95 | + |
| 96 | + // for (feat <- src.features) { |
| 97 | + // val intersections = |
| 98 | + // src.filter(filter.Filter.intersects(feat.geometry)) |
| 99 | + // dest ++= |
| 100 | + // intersections.filter(_.id > feat.id).map { corner => |
| 101 | + // feature.Feature( |
| 102 | + // "geom" -> (feat.geometry intersection corner.geometry), |
| 103 | + // (joinField + "Left") -> feat.get[Any](joinField), |
| 104 | + // (joinField + "Right") -> corner.get[Any](joinField) |
| 105 | + // ) |
| 106 | + // } |
| 107 | + // } |
| 108 | + |
| 109 | + // println("Found %d intersections".format(dest.count)) |
| 110 | + // } |
| 111 | + |
| 112 | + // def rewrite(schema: feature.Schema, fieldName: String): feature.Schema = |
| 113 | + // feature.Schema( |
| 114 | + // schema.name + "_intersections", |
| 115 | + // feature.Field( |
| 116 | + // "geom", |
| 117 | + // classOf[com.vividsolutions.jts.geom.Geometry], |
| 118 | + // schema.geometry.projection |
| 119 | + // ), |
| 120 | + // feature.Field(fieldName + "Left", classOf[String]), |
| 121 | + // feature.Field(fieldName + "Right", classOf[String]) |
| 122 | + // ) |
| 123 | + |
| 124 | + // def main(args: Array[String]) = { |
| 125 | + // if (args.length == 0) { |
| 126 | + // println("You need to provide the path to a shapefile as an argument to this example.") |
| 127 | + // } else { |
| 128 | + // val src = layer.Shapefile(args(0)) |
| 129 | + // val joinField = |
| 130 | + // src.schema.fields.find { _.gtBinding == classOf[String] } match { |
| 131 | + // case Some(f) => f.name |
| 132 | + // case None => "id" |
| 133 | + // } |
| 134 | + // val dest = src.workspace.create(rewrite(src.schema, joinField)) |
| 135 | + // process(src, dest, joinField) |
| 136 | + // } |
| 137 | + // } |
51 | 138 | }
|
0 commit comments