Documentation was a bit confusing

I'm not a GIS guru, I'm just posting this to provide some feedback on the documentation. (the struggles of a beginner)

This is just for reference and in case anyone else stubs their toes. Or if I forget how to do this and need to look back later :) Feel free to close this issue

For background:

I took the world shorelines from

https://www.ngdc.noaa.gov/mgg/shorelines/

I converted it to a GeoJSON using ogr2ogr

gr2ogr -f GeoJSON shoreline.json GSHHS_c_L1.shp

I then read it into Clojure with

(-> "/path-to-my-file/shoreline.json"
    slurp
    read-geojson)

So far so good :)

I then want to cut down these shoreline polygons to a certain lat-long region of interest. (I will later convert this to an SVG). So I try to create a polygon and then intersect it with my shoreline

And this is where things get hairy. I try to create a polygon

http://factual.github.io/geo/3.0.1/geo.jts.html#var-polygon

the docs read

(polygon shell)
(polygon shell holes)

Given a LinearRing shell, and a list of LinearRing holes, generates a
polygon.

Okay, so to make a polygon I need to first make a LinearRing. I looke up the docs for that: http://factual.github.io/geo/3.0.1/geo.jts.html#var-linear-ring

linear-ring
(linear-ring coordinates)
(linear-ring coordinates srid)

Given a list of Coordinates, creates a LinearRing. Allows an optional SRID argument at end.

Okay so to make a linear ring I need to make a Coordinates ... I look up that: http://factual.github.io/geo/3.0.1/geo.jts.html#var-coordinates

coordinates
(coordinates geom)

Get a sequence of Coordinates from a Geometry

Now i'm way down the rabbit hole, but there is no geometry

There are two similar things

geometries
(geometries c)

Given a GeometryCollection, generate a sequence of Geometries

and

geometry-collection
(geometry-collection geometries)

Given a list of Geometries, generates a GeometryCollection.

But as you can see they reference each other in their arguments

So it turns out that when it says (linear-ring coordinates) the docs read "Given a list of Coordinates ..." what this expects is an actual Clojure seq with coordinates.

First I tried to do

(linear-ring [ (coordinate 0 0) (coordinate 0 75) (coordinate 75 75) (coordinate 75 0) ])

And I got an error. Turns out you need to make sure the last point is equal to the first

So this works

(linear-ring [ (coordinate 0 0) (coordinate 0 50) (coordinate 50 50) (coordinate 50 0)  (coordinate 0 0) ]

Now that I've made a polygon I can intersect it

Putting it all together I got this:

  (let [shoreline   (-> "/my-path-to-shoreline-file/shoreline.json"
            slurp
            read-geojson)
    area-of-interest (polygon (linear-ring [(coordinate 0 0)
                        (coordinate 0 50)
                        (coordinate 50 50)
                        (coordinate 50 0)
                        (coordinate 0 0) ]))]
    (->> shoreline
     (map (fn [shoreline-poly]
        (update shoreline-poly
            :geometry
            #(intersection %
                       area-of-interest))))
      (filter #(not (.isEmpty (:geometry %))))
     to-geojson-feature-collection
     (spit "trimmed-shorelines.json")))

And to my surprise trimmed-shorelines.json opens up in QGIS. And it even looks correct!

(Please let me know if I'm not doing this in the way the library expects me to :) )

The only caveat is that at certain "clip" sizes I'd get two layers in the GeoJSON for some reason. Most would be Polygons and once in a while I'd get GeometryCollection. My best guess is one of the intersections split a (concave) polygon into multiple smaller polygons and this gets crammed into a single collection.

If my guess is right then I imagine with some fiddling you could probably duplicate the Shapelike that has the GeometryCollection and copy out each Polygon to flatten the whole thing - so that you only have Shapelike's with polygons

I haven't done the Polygon to SVG conversion yet. I'm not super clear on how the library expects me to extract coordinates. The only way I see for now is to dip into JTS. It's a tad ugly but it seems to work.

  (let [shoreline   (-> "/home/geokon/Downloads/coast/GSHHS_shp/c/shoreline.json"
            slurp
            read-geojson)
    area-of-interest (polygon (linear-ring [(coordinate 0 0)
                        (coordinate 0 50)
                        (coordinate 50 50)
                        (coordinate 50 0)
                        (coordinate 0 0) ]))]
      (->> shoreline
       (map (fn [shoreline-poly]
          (update shoreline-poly
              :geometry
              #(intersection %
                     area-of-interest))))
       (filter #(not (.isEmpty %)))
     (map :geometry)
     (map coordinates)
     (map coord-array) 
     (map (fn [poly]
        (map #(vector (.getX %)
                  (.getY %))
             poly)))

Please let me know if there is something more kosher

Once I get a basic plotting/display pipeline going I'll post an update. Thanks for the library guys! It's going way smoother than I anticipated (my only hangup really was how to read in shapefiles)

Okay, with a little futzing around I managed to massage the output into an SVG

Here is the result: https://geokon-gh.github.io/map-demo.html

(scroll to bottom to see the picture :) )

The unrolling of the geometries into polygons isn't exactly robust - but it was good enough for the moment

EDIT: As a sanity check I changed it from drawing svg/polygon to svg/line-strip and it still looks good

@geokon-gh thanks for documenting these questions you came across. I think the key issue you have hit upon is what you summarized here:

So it turns out that when it says (linear-ring coordinates) the docs read "Given a list of Coordinates ..." what this expects is an actual Clojure seq with coordinates.

There are different APIs in the library depending on whether you're working with native JTS types (like org.locationtech.blah.Polygon) or the Clojure-ish representations of these structures which mirror what is used in things like WKT or GeoJSON (so a point is [1 1] and a linestring is [[1 1] [0 0]] etc)

We should do a better job of documenting these, probably the best way would be by adding more examples into the docstrings.

1 thing that might be useful for you is the variety of blah-wkt functions, which are designed to turn clojure-ish nested vector structures into the appropriate JTS type, e.g.:

linestring-wkt
(linestring-wkt coordinates)
(linestring-wkt coordinates srid)

Makes a LineString from a WKT-style data structure: a flat sequence of
coordinate pairs, e.g. [0 0, 1 0, 0 2, 0 0]. Allows an optional SRID argument at end.

So I think you might be able to do something like

;; hope I've got the bracket balancing right here, but you get the idea
(polygon-wkt [[[0 50] [50 50] [50 0] [0 0] [0 50]]])

Sorry for the slow reply. I just wanted to take the time to look at this properly.

Thank you for showing me how to do it more correctly. I actually didn't realize there is a while Clojure abstraction layer going on. I guess what threw me off was the function signature. I only knew that WKT is some other format I'm not interested in so I didn't even bother looking into it haha :smile: - but what you're showing me makes a ton of sense b/c JSON is like an ugly EDN. So you'd want to represent all the Geometry as Clojure datastructures when possible.

I'm still trying to wrap my head around some of the design decisions. I guess my question would be, why is the JTS layer exposed so much? When i do read-geojson why am I getting opaque Java/JTS objects in my map for the Polygons? Wouldn't it make more sense to get back a straight JSON->EDN conversion with Clojure maps/vectors. In-effect you'd have a GeoEDN

Then those geometry vectors would directly correspond to what you are passing into polygon-wkt in your example (well the function would be redundant I guess).

At the next step I do this intersection between two geometries. It's a bit weird how I need to dig into the guts and intersect the :geometry values from each object - but maybe this is unavoidable. I'm guessing that the other fields ( :properties ?) don't have a logical intersection?

I also get how this is really calling out to some JTS functions, so having JTS objects as input does seemingly simplify things. Maybe this is the answer to my previous question :) But I think you could just do the conversion from EDN to JTS objects under the hood. It'd involve a bit more machinery, but if you memoized the conversion then it'd just happen once for each polygon vector. You could even run the memoized conversion during read-in.

The only catch would be that each time the user tweaked the EDN it'd have to be reconverted (while if they somehow modify the JTS object directly then it doesn't need to worry about that), but I don't think that's too unreasonable..

Edit: Anyway, those are maybe some disjoint thoughts. I don't mean to be critical - just trying to grasp the design. This library is great and helping me a lot :) You can do something vaguely similar with SVG polygons:

https://github.com/thi-ng/geom/blob/master/geom-svg/src/core.org#example-use-case

https://github.com/thi-ng/geom/blob/85783a1f87670c5821473298fa0b491bd40c3028/org/src/types/polygon.org

It uses the "Sutherland–Hodgman algorithm" - however I honestly haven't tested it myself - and I'd assume it's not as performant of course

And SVG polygons are annoyingly slightly different from GeoJSON ones (you can't cut holes in the same way). And furthermore, the ones in geom don't seem to make hole at all (I'm probably misunderstanding the design of geom though). Its not an issue in my scenario, but it's a limitation for a universal one-to-one correspondence

I'm still trying to wrap my head around some of the design decisions. I guess my question would be, why is the JTS layer exposed so much? When i do read-geojson why am I getting opaque Java/JTS objects in my map for the Polygons? Wouldn't it make more sense to get back a straight JSON->EDN conversion with Clojure maps/vectors. In-effect you'd have a GeoEDN

Well the main reason is for efficiency...JTS is where the actual geometry operations are all implemented (and optimized) and the GeoJSON-in-EDN formats are really just there for dealing with serialization, and for convenience like if you quickly need to just construct a polygon by hand from the repl or something.

At the next step I do this intersection between two geometries. It's a bit weird how I need to dig into the guts and intersect the :geometry values from each object - but maybe this is unavoidable. I'm guessing that the other fields ( :properties ?) don't have a logical intersection?

GeoJSON can be a little confusing and this structure with :properties and :geometry is designed to encode a "Feature" that has additional metadata associated with it.

If you don't have that case (i.e. you only have geometry) you might look at some of the other functions in the io namespace -- there should be one to read only a geometry.

I'm guessing that the other fields ( :properties ?) don't have a logical intersection?

That sounds right to me; I think you'd have to implement your own logic for how you want to handle these. Depending on your use case you might want to just keep one or the other, or merge them, etc.

Thanks for clarifying and helping me grok the logic a bit :)

Going to close this as I think we hopefully sorted out these issues in the previous discussion. Feel free to reopen (or file a new issue) if you have further questions.

Factual / geo

Documentation was a bit confusing #71