Geometry::get_point_vec() returns empty Vec

[metasim]

I'm unable to access get the points of any GDAL Geometry instance I create. The following will fail at the last test:

#[test]
fn test_get_points() -> std::result::Result<(), GdalError> {
    let sample = "/vsicurl/https://gist.githubusercontent.com/metasim/15fb05bbe76c1d7e5583033f2b6757f7/raw/fac8b385b733f50f0b7a1466db03e7918084b4a1/test.geojson";
    let ds = Dataset::open(sample)?;
    let layer = ds.layer(0)?;
    let feature = layer.feature(0).unwrap();
    let geom = feature.geometry();

    assert_eq!(geom.geometry_count(), 1);
    assert!(geom.area() > 0.);
    assert_eq!(geom.geometry_type(), gdal_sys::OGRwkbGeometryType::wkbPolygon);
    assert!(geom.json()?.contains("Polygon"));
    let points = geom.get_point_vec();
    assert!(!points.is_empty());

    Ok(())
}

[metasim]

An even simpler test case:

#[test]
fn test_get_points() -> std::result::Result<(), GdalError> {
    let geom = Geometry::bbox(0., 0., 1., 1.)?;
    assert_eq!(geom.geometry_count(), 1);
    assert!(geom.area() > 0.);
    assert_eq!(geom.geometry_type(), gdal_sys::OGRwkbGeometryType::wkbPolygon);
    assert!(geom.json()?.contains("Polygon"));

    let points = geom.get_point_vec();
    assert!(!points.is_empty());

    Ok(())
}

[metasim]

Interestingly, running test_geom_accessors passes, but pasting the above code into the same file and running it still fails.

gdal/src/vector/vector_tests/mod.rs

Lines 485 to 509 in 492e0c8

	fn test_geom_accessors() {
	with_feature("roads.geojson", 236194095, |feature| {
	let geom = feature.geometry();
	assert_eq!(geom.geometry_type(), OGRwkbGeometryType::wkbLineString);
	let coords = geom.get_point_vec();
	assert_eq!(
	coords,
	[
	(26.1019276, 44.4302748, 0.0),
	(26.1019382, 44.4303191, 0.0),
	(26.1020002, 44.4304202, 0.0)
	]
	);
	assert_eq!(geom.geometry_count(), 0);
	let geom = feature.geometry_by_index(0).unwrap();
	assert_eq!(geom.geometry_type(), OGRwkbGeometryType::wkbLineString);
	assert!(feature.geometry_by_index(1).is_err());
	let geom = feature.geometry_by_name("");
	assert!(geom.is_ok());
	let geom = feature.geometry_by_name("").unwrap();
	assert_eq!(geom.geometry_type(), OGRwkbGeometryType::wkbLineString);
	assert!(feature.geometry_by_name("FOO").is_err());
	});
	}

[metasim]

Would love to know if I'm doing something wrong, or if there's a workaround (e.g. call methods in certain order).

[lnicola]

The unit test uses a line string. Polygons don't contain points, but rings, so you need to go through those. I can't give an example right now.

[metasim]

I don't see any methods to get rings. Only points.

I took a look at the implementation of impl TryFrom<Geometry> for geo_types::Geometry<f64> and saw that it uses the unsafe get_unowned_geometry method to dig down into the lower structures. I can get what I want with the following call:

unsafe { geom.get_unowned_geometry(0) }.get_point_vec()

However, this doesn't really seem like the best user-facing method.

I'd like to submit a PR for something more friendlly, but don't want to waste time on something that isn't inline with the maintainers' thinking.

From my standpoint, modifying the Geometry::get_point_vec() method to got into the unowned geometries automatically seems like better ergonomics. Does that approach pass your sniff test, or should I take a different approach (like, have methods to expose rings)?

[lnicola]

From my standpoint, modifying the Geometry::get_point_vec() method to got into the unowned geometries automatically seems like better ergonomics. Does that approach pass your sniff test, or should I take a different approach (like, have methods to expose rings)?

I don't think it makes sense to return points from all the rings together. I think the first ring is the exterior of the polygon, while the others are "holes".

We should try match the GDAL API:

>>> ring = ogr.Geometry(ogr.wkbLinearRing)
>>> ring.AddPoint(1179091.1646903288, 712782.8838459781)
>>> ring.AddPoint(1161053.0218226474, 667456.2684348812)
>>> ring.AddPoint(1214704.933941905, 641092.8288590391)
>>> ring.AddPoint(1228580.428455506, 682719.3123998424)
>>> ring.AddPoint(1218405.0658121984, 721108.1805541387)
>>> ring.AddPoint(1179091.1646903288, 712782.8838459781)
>>> poly = ogr.Geometry(ogr.wkbPolygon)
>>> poly.AddGeometry(ring)
0
>>> poly.GetPointCount()
0
>>> poly.GetPoints()
ERROR 6: Incompatible geometry for operation
>>> ring = poly.GetGeometryRef(0)
>>> ring.GetGeometryName()
'LINEARRING'
>>> ring.GetPoints()
[(1179091.1646903288, 712782.8838459781, 0.0), (1161053.0218226474, 667456.2684348812, 0.0), (1214704.933941905, 641092.8288590391, 0.0), (1228580.428455506, 682719.3123998424, 0.0), (1218405.0658121984, 721108.1805541387, 0.0), (1179091.1646903288, 712782.8838459781, 0.0)]

Notice how:

• GetPoints is fallible
• you can get a reference to an owned geometry using GetGeometryRef
• our API isn't really idiomatic, hence the unsafe get_unowned_geometry and the don't keep this object for long comment in its implementation; the return type is logically a &'self Geometry, since the ring is owned by the polygon.

[metasim]

We should try match the GDAL API:

Very helpful to know where to go for a "reference implementation". Thanks! :-)

[rmanoka]

Also, see geo_to_gdal.rs and gdal_to_geo.rs. Those two files were an awesome reference to accessing inner structures using our bindings.

[metasim]

What is the definition of "(un)owned" in get_unowned_geometry? Is a geometry "owned" when:

a. The lifecycle is managed by Rust; if the binding leaves scope, then Drop is invoked. Or...
b. The lifecycle is managed by gdal; if the binding leaves scope, we assume gdal will deallocate when appropriate.

My real question is how to think about returning a reference to something that is logically bound to 'self, but across the FFI boundary. I'm new to unsafe Rust, so the best I've come up with so far is this (which SIGSEGVs, and shows my ignorance of the pointer semantics):

    /// Get a reference to the contained geometry at `index`.
    pub fn get_geometry<'a>(&'a self, index: usize) -> Result<&'a Geometry> {
        let ring = unsafe { self.get_unowned_geometry(index) };
        let ring_ref: &'a Geometry  = unsafe { &*(&ring as * const Geometry) };
        Ok(ring_ref)
    }

(I also tried using std::mem::transmute with same result.)

I'm also wondering if I should be considering a GeometryRef type (like the Python example), using PhantomData to track lifetime?

[metasim]

To further refine my questions, I'm interpreting this:

gdal/src/vector/geometry.rs

Lines 306 to 308 in 8cb6e92

	pub unsafe fn get_unowned_geometry(&self, n: usize) -> Geometry {
	// get the n-th sub-geometry as a non-owned Geometry; don't keep this
	// object for long.

and this:

our API isn't really idiomatic, hence the unsafe get_unowned_geometry and the don't keep this object for long comment in its implementation;

To mean the naïve implementation such as this will eventually break:

    pub fn get_geometry(&self, index: usize) -> Result<Geometry> {
        let ring = unsafe { self.get_unowned_geometry(index) };
        Ok(ring)
    }

With that assumption, I've given up on figuring out how to return Result<&Geometry> and am taking the approach of using a wrapper type.

The implementation below works. However, I'm not confident enough in my 'lifetime skills to know if I'm just getting lucky with the test, or if this is indeed a proper way of returning a Geometry that is tied to the lifetime of self.

pub struct GeometryRef<'a> {
    geom: Geometry,
    _lifetime: PhantomData<&'a ()>
}

impl <'a> Deref for GeometryRef<'a> {
    type Target = Geometry;

    fn deref(&self) -> &Self::Target {
        &self.geom
    }
}

...
    pub fn get_geometry(&self, index: usize) -> Result<GeometryRef> {
        let ring = unsafe { self.get_unowned_geometry(index) };
        let gref = GeometryRef { geom: ring, _lifetime: PhantomData::default() };
        Ok(gref)
    }

...
    // This (mostly) works using the `GeometryRef` approach.
    #[test]
    fn test_ring_points() {
        let mut ring = Geometry::empty(wkbLinearRing).unwrap();
        ring.add_point_2d((1179091.1646903288, 712782.8838459781));
        ring.add_point_2d((1161053.0218226474, 667456.2684348812));
        ring.add_point_2d((1214704.933941905, 641092.8288590391));
        ring.add_point_2d((1228580.428455506, 682719.3123998424));
        ring.add_point_2d((1218405.0658121984, 721108.1805541387));
        ring.add_point_2d((1179091.1646903288, 712782.8838459781));
        assert!(!ring.is_empty());
        assert_eq!(ring.get_point_vec().len(), 6);
        let mut poly = Geometry::empty(wkbPolygon).unwrap();
        poly.add_geometry(ring.to_owned()).unwrap();
        // Points are in ring, not containing geometry.
        // TODO: In Python SWIG bindings, `GetPoints` is fallible.
        assert!(poly.get_point_vec().is_empty());
        assert_eq!(poly.geometry_count(), 1);
        let ring_out = poly.get_geometry(0).unwrap();
        println!("{}", ring_out.wkt().unwrap());
        // TODO: This fails. Returns `wkbLineString`, but `wkb()` shows it to be a `LINEARRING`
        //assert_eq!(ring_out.geometry_type(), wkbLinearRing);
        assert!(!&ring_out.is_empty());
        assert_eq!(ring.get_point_vec(), ring_out.get_point_vec());
    }

[metasim]

@lnicola @rmanoka Wondering if either of you have advice to share on the appropriate route to take with the above. I think the GeometryRef approach is working.

[lnicola]

Sorry, I've been less around lately. It's not really clear to me why the original approach didn't work. GeometryRef might work -- it takes the lifetime of self, right?

[metasim]

Sorry, I've been less around lately.

Well, sorry to bother you!. I'm excited this project exists at all, and want to help out when I run into friction. That said, I'm still struggling to understand lifetimes when FFI is concerned, and may require more hand-holding than you're able to provide. IOW, helping out may be above my current Rust skillset.

it takes the lifetime of self, right?

Yes.

[pka]

Reading Polygon points like this should work:

let ring_count = geom.geometry_count();
for i in 0..ring_count {
    let ring = unsafe { geom.get_unowned_geometry(i) };
    let length = unsafe { gdal_sys::OGR_G_GetPointCount(ring.c_geometry()) } as usize;
    for i in 0..length {
        let (x, y, z) = ring.get_point(i as i32);
    }
}

[metasim]

@pka Thanks, I've seen that code. I'm actually want the whole Geomtery type, not just the points.