ewestern / geos

This is a Haskell binding to Geos, the open-source geometry library
MIT License
13 stars 9 forks source link

Null pointer when trying to pull out CoordSeq of a Geometry Polygon #35

Closed durgaBahadur closed 3 years ago

durgaBahadur commented 3 years ago

Hi there!

First off, thanks for this great library!

I am having a hard time trying to extract the coordinates out of a 'Geometry Polygon'. Unfortunately I get a null pointer exception. Please see the minimal ghci-demo I've created to illustrate the problem:

Prelude> :l GetCoordSeqDemo.hs
[1 of 1] Compiling GetCoordSeqDemo  ( GetCoordSeqDemo.hs, interpreted )
Ok, one module loaded.
*GetCoordSeqDemo> testCoordSeqGetter $ getPgRawGeom  createPg
*** Exception: Encountered null pointer at: getCoordinateSequence
CallStack (from HasCallStack):
  error, called at src/Data/Geometry/Geos/Raw/Base.hs:68:14 in geos-0.4.1-ABmO5CcGEiV20rMbTpIAA0:Data.Geometry.Geos.Raw.Base

And this would be the content of GetCoordSeqDemo.hs

{-# LANGUAGE DataKinds           #-}
{-# LANGUAGE FlexibleInstances   #-}
{-# LANGUAGE GADTs               #-}
{-# LANGUAGE OverloadedStrings   #-}
{-# LANGUAGE ScopedTypeVariables #-}

module GetCoordSeqDemo where

import qualified Data.Vector                     as V

import           Data.Geometry.Geos.Geometry
import qualified Data.Geometry.Geos.Raw.Base     as RB
import qualified Data.Geometry.Geos.Raw.CoordSeq as RC
import qualified Data.Geometry.Geos.Raw.Geometry as R

-- stolen from Helpers.hs
makeLinearRing :: [(Double, Double)] -> Either GeometryConstructionError LinearRing
makeLinearRing = linearRing . V.map (uncurry coordinate2) . V.fromList

makePolygon :: [[(Double, Double)]] -> Either GeometryConstructionError Polygon
makePolygon ls = do
  lRings <- traverse makeLinearRing ls
  polygon $ V.fromList lRings

makePolygonGeo :: [[(Double, Double)]] -> Either GeometryConstructionError (Geometry Polygon)
makePolygonGeo cs = (\pg -> PolygonGeometry pg Nothing) <$> makePolygon cs

fromRight :: Either a b -> b
fromRight (Right v) = v
fromRight _         = error "Shouldn't happen."

-- demo starts here
createPg :: Geometry Polygon
createPg = (fromRight . makePolygonGeo) [[(1.1,1.2),(22.1,22.2),(33.1,33.2),(1.1,1.2)]]

getPgRawGeom :: Geometry Polygon -> R.GeomConst
getPgRawGeom geometryPg = RB.runGeos $ do
  geometryP' :: R.GeomConst <- convertGeometryToRaw geometryPg
  return geometryP'

testCoordSeqGetter :: R.Geometry a => a -> RC.CoordSeq
testCoordSeqGetter geom = RB.runGeos $ do
  R.getCoordinateSequence geom

Since I am relatively new to Haskell I am not sure if I am doing sth. terribly wrong here or if I found a bug. Any help would be appreciated!

ewestern commented 3 years ago

So, first thing to note is that you're using the Raw interface. This is fine, but it's mostly only exposed for people who prefer to use a low-level interface for performance reasons. If you don't have a strong reason for using it, then I would suggest using modules without Raw in the name.

The reason this is happening is because when you use the Raw interface, the geometry objects you are using are essentially thin wrappers around c structs, along with some logic to specify when to garbage collect those objects. A call to runGeos is meant to orchestrate the memory management among all these different objects; as such, you aren't meant to return any Raw objects from a call to runGeos, because, essentially, the memory representing that object will already be freed once runGeos completes.

This should all be handled for you if you use, for example, Geometry Polygon, and just run whatever computation you want to run on that value, rather than converting it to a Raw object.

durgaBahadur commented 3 years ago

Cool, thanks for the explanation! The reason I found myself digging in the raw is because I need to get my hands on the bare x/y coordinates of my polygons. At top level I couldn't find anything that would do that, so I hoped something like this might give me the coordinate pairs:

pgToCoordList :: Geometry Polygon -> [(Double,Double)]
pgToCoordList geomPg = RB.runGeos $ do
  geomPg' :: R.GeomConst <- convertGeometryToRaw geomPg
  cs <- R.getCoordinateSequence geomPg'  -- exception
  dim <- RC.getCoordinateSequenceSize cs
  mapM (extractDouble cs) [0 .. (dim - 1)]
  where extractDouble cs idx = do
          x <- RC.getCoordinateSequenceX cs idx
          y <- RC.getCoordinateSequenceY cs idx
          return (x,y)

With these coordinates I plan to phrase a query for an api that exclusively accepts Esri geometry syntax as parameters.

Thirst thought was to do pattern matching up to the points, but therefore I would need all the data constructors to be exported which isn't the case. Possibly there is an easier method I have overlooked?

ewestern commented 3 years ago

I was about to write out a simple function for that, but I realized that one of the constructors I need to do that isn't exposed, which is probably why you went looking in the Raw interface. My fault. Will expose a constructor and record name, which should make it much more straightforward.

durgaBahadur commented 3 years ago

Nice! Actually, I could do that and send a PR. Would this be helpful?

ewestern commented 3 years ago

Please do. It might be a few days before I can get around to it.

durgaBahadur commented 3 years ago

Great, I won't be able to continue work this week. But next or latest after next week will be possible, and I am glad to contribute!