Sg/area dist debug

[skygering]

I realized that I missed a few geometry types that were coming up in my simulations.

I also realized that my code failed with empty geometries as it was trying to access the first point. This then led me to the discovery of an issue with the implementation of isempty at least with LibGEOS. So now I have to check isempty and find if there is a non-zero number of points to determine the type of the points within a geometry that returns an area of 0.

import LibGEOS as LG
import GeoInterface as GI

p = LG.readgeom("POLYGON EMPTY")
GI.isempty(p)  # false
GI.npoint(p) # 0

This seems like too much work just to return 0. I am wondering if it would be bad to just return an Int 0 for everything but non-empty polygons / multipolygon / collections.

This isn't the best type-wise, but there is no guarantee that we can even determine the type.

Thoughts?

[rafaqz]

I think we need to either pick a fixed floating point output type (e.g. Float64), or detect it once very early and push that through. Then all those zeros will be type stable and we cont need to do the zero(typeof(x(...

[rafaqz]

Could we make area always return Float64 for type stability?

Another option is to check the type once at the start, then push it through as the first argument T to an _area(T, trait, geom) method that returns e.g. zero(T). We can let users specify it too if they need to

[rafaqz]

Same here

[skygering]

I wanted to preserve the flexibility so that if I run a simulation with Float32 or Float64 polygons, the area returns the correct type. But I can see why that is tricky, especially with the empty geometries for which we would have to default to something. Do you think it is better just for it to return Float64 so that the user can plan?

[rafaqz]

Maybe default to Float64 but let the user pass in the type?

[skygering]

Okay, I updated the type stuff! They all have optional parameters now. I also just simplified some stuff.

Something I am noting is that I am doing GI.getpoint(geom, i) a lot in order to deal with "closed" curves. This is my way of making it so that it doesn't matter if a Geometry library repeats the last coordinates in linear rings / polygons, we still repeat the last point for our algorithms. I don't have an idea for a better way of doing this though.

As discussed in #32, this can be problematic if a library does a poor job of implementing GeoInterface. I kinda feel like that is the library's fault though, not ours. For example, using tuples to convert an input geometry into a GeoInterface geometry feels like cheating, because we aren't really dealing with an arbitrary polygon, we're converting a given polygon to a specific type...

However, if you think this is still slower than GI.getpoint(geom) and using the iterator regardless, then we can think up a new way of iterating. Maybe that can be a separate PR though.

[rafaqz]

Looks like it will be much faster :)

Lets not worry about getgeom for now. But I do worry that getting individual points from LibGeos or ArchGDAL will never be fast (and geointerface should actually caution against operating at the point level too much)