Re: [gdal-dev] GDAL, ESRI Shapefiles and nested polygons
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
> To: Graeme Wilkie
> Cc: "gdal-dev@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
> >
> > To: gdal-dev@lists.osgeo.org
> > Cc: Graeme Wilkie
> > 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 fea
Re: [gdal-dev] GDAL, ESRI Shapefiles and nested polygons
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
To: Graeme Wilkie
Cc: "gdal-dev@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
> To: gdal-dev@lists.osgeo.org
> Cc: Graeme Wilkie
> 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/___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
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
> To: gdal-dev@lists.osgeo.org
> Cc: Graeme Wilkie
> 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
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
Re: [gdal-dev] GDAL, ESRI Shapefiles and nested polygons
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
If anybody has a small piece of example could it would be a great help.
Regards
Graeme
From: Even Rouault
To: gdal-dev@lists.osgeo.org
Cc: Graeme Wilkie
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/___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
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 ___ gdal-dev mailing list gdal-dev@lists.osgeo.org http://lists.osgeo.org/mailman/listinfo/gdal-dev
[gdal-dev] GDAL, ESRI Shapefiles and nested polygons
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.
If I use the lower level routines SHPOpen, SHPReadObject, SHPClose I can get
access to the data (after a fashion) but the requirement is to use the
higher level API.
I'm using GDAL-1.6.2 (yes I know I should be using a much more up-to-date
version but the project will not sanction it). I'm developing on a Linux
(Redhat 5.9) system in GCC.
I'd appreciate it if you could point out what I'm doing wrong (I suspect it
is my error).
OGRDataSourceH hDS;
OGRFeatureH hFeature;
OGRLayerH;
OGRRegisterAll();
hDS = OGROpen(filename, FALSE, NULL);
hLayer = OGR_DS_GetLayer(hDS, 0); // only interested in layer 0
while ((hFeature = OGR_L_GetNextFeature(hLayer)) != NULL)
(
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
do
{
OGRGeometryH hRing = OGR_G_GetGeometryRef(hGeometry, index);
Int num_vertices = OGR_G_GetPointCount(hRing); // returns 0
//.
//. Do something with the data
//.
Index++;
}
// Destroy the hfeature here
)
___
gdal-dev mailing list
gdal-dev@lists.osgeo.org
http://lists.osgeo.org/mailman/listinfo/gdal-dev
