[gdal-dev] GDAL, ESRI Shapefiles and nested polygons

[Even Rouault]

Le vendredi 08 août 2014 16:54:04, Graeme Wilkie a écrit :
> Hi,
>  
> I don't have a WKT format as such (sorry but I'm learning about all this as
> I go), all I do have is a simple Python script that creates the
> shapefiles. The idea of having T3 inside T2 with both inside /t1 is
> correct in this case but not necessarily always the case. T2 & T3 could
> both be inside T1 side by side. 
> The content of the python script file is as follows (the polygons are T1,
> T2, T3 in order): 
> import shapefile
> outfile = shapefile.Writer(shapefile.POLYGON)
> outfile.poly(parts=[\
> [[easting,northing],[easting,northing],[easting,northing],[easting,northing
> ],\
> [easting,northing],[easting,northing],[easting,northing],[easting,northing
> ]],\
> [[easting,northing],[easting,northing],[easting,northing],[easting,northin
> g]]\ ])
> outfile('Location','C','20')
> outfile.record('One')
> outfile.save('nestedpolygon')
>  
> I've replaced to values with northing and easting to save having to enter
> the figures (I can't include the file). 
> I'm sure my mistake is obvious but as I'm still pretty much a beginner I
> can't see it :( 

What you must realize is that the models of shapefiles and OGR geometry (the 
latter following simple feature specifications) are different. shapefiles have a 
"polygon" type that might be in OGR term either a polygon or a multipolygon. 
So you should be able to handle both cases.
A OGR polygon is a collection of rings, the first one being the outer ring and 
the following ones the inner rings
A OGR multipolygon is a collection of polygon.
When reading a shapefile, OGR does some topological analysis to figure out the 
spatial relationship between rings and determine if they must be outer or 
inner rings.

> All suggestions are welcome.
> 
> ________________________________
>  From: Even Rouault <even.rouault at spatialys.com>
> To: Graeme Wilkie <graeme.wilkie at btinternet.com>
> Cc: "gdal-dev at lists.osgeo.org" <gdal-dev at lists.osgeo.org>
> Sent: Friday, August 8, 2014 2:27 PM
> Subject: Re: [gdal-dev] GDAL, ESRI Shapefiles and nested polygons
> 
> Le vendredi 08 août 2014 15:19:38, Graeme Wilkie a écrit :
> > Hi,
> >
> > 
> >
> > You were correct, it is returning multipolygon (9sorry, I should have
> > included that information).
> > I've made the change you suggested so the code looks like
> >
> >  Int index = 0;
> >
> > OGRGeometryH hGeometry = OGR_F_GetGeometryRef(hFeature);
> > int shape_type = wkbFlatten(ORG_G_GetGeometry(hGeometry);
> > int num_geom = OGR_G_GetGeometryCount(hGeometry) - 1; // Returns a valid
> > result int innerIndex = 0;
> >
> > 
> >
> > do
> > {
> >
> >                 OGRGeometryH hRing = OGR_G_GetGeometryRef(hGeometry,
> >
> > index); OGRGeometryH tmpGeo = OGR_G_GetGeometryRef(hRing, innerIndex);
> > Int num_vertices = OGR_G_GetPointCount(tmpGeo);
> >
> >                 //.
> >
> > //. Do something with the data
> > //.
> > Index++;
> > }
> > // Destroy the hfeature here
> > )
> >
> > 
> >
> > If innerIndex is 0 then I can read the data for that polygon. If
> > innerIndex in 1 then tmpGeo is NULL even though there are 2 inner
> > polygons in the test data.
> > Have I do something stupid ?
> >
> > 
> >
> > The test data is a simple square with 2 squares inside it created using
> > shapefile.py
> 
> You still didn't include the WKT... Well I will assume that you have square
> T3 included inside square T2 included inside square T1.
> 
> In that case the structure of the WKT should be ((T1,T2),T3)
> i.e. a multipolygon with 2 polygons : (T1,T2) and T3
> T1 is an outer ring
> T2 an inner ring of T1
> T3 an outer ring
> That's the way polygons are modelised in SFS specification.
> 
> > If anybody has a small piece of example could it would be a great help.
> >
> > 
> >
> > Regards
> >
> > 
> >
> > Graeme
> > 
> > ________________________________
> >
> >  From: Even Rouault <even.rouault at spatialys.com>
> >
> > To: gdal-dev at lists.osgeo.org
> > Cc: Graeme Wilkie <graeme.wilkie at btinternet.com>
> > Sent: Thursday, August 7, 2014 9:26 PM
> > Subject: Re: [gdal-dev] GDAL, ESRI Shapefiles and nested polygons
> > 
> > Le jeudi 07 août 2014 22:20:34, Graeme Wilkie a écrit :
> > > Hi,
> > > 
> > > I'm trying to use GDAL to read ESRI shapefiles. I've got it working for
> > > points, lines, polygons, polylines and multipoint but I'm having
> > > problems with polygons that contain polygons in the same feature.
> > > 
> > > 
> > > 
> > > The following is a cut down code sample of the code I'm using. The
> > > problem I have is that I can read the number of polygons in the feature
> > > (i.e. 3 polygons, nested) using OGR_G_GetGemoetryCount. When I use
> > > OGR_R_GetGeometryRef to get the handle/pointer for each of the polygons
> > > in turn I get what appears to be a valid handle but when I use it to
> > > get the number of points/vertices in the polygon it always returns 0.
> > 
> > Graeme,
> > 
> > You didn't display the WKT that corresponds to the geometry, so I will
> > just make the guess that the geometries returned is a multipolygon and
> > not a polygon.
> > 
> > So you might need one more level of OGR_G_GetGeometryCount and
> > OGR_G_GetGeometryRef to go to the ring level.
> > 
> > Best regards,
> > 
> > Even

-- 
Spatialys - Geospatial professional services
http://www.spatialys.com