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

[Graeme Wilkie]

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,northing]]\
])
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 :(
 
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/
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.osgeo.org/pipermail/gdal-dev/attachments/20140808/151559ee/attachment.html>