 |
||||
|
||||
[Shapelib] Identify polygon by XY coordinate?shapelib@lists.maptools.org shapelib@lists.maptools.orgWed, 17 Mar 2004 10:37:35 -0500
|
||||
This is a multi-part message in MIME format.
--------------030000040708090904020500
Content-Type: text/plain; charset=us-ascii; format=flowed
Content-Transfer-Encoding: 7bit
shapelib-admin@lists.maptools.org wrote:
> List,
> I need a function that uses (X,Y) coordinates to identify the
> surrounding polygon in a shapefile, and then retrieve the associated
> attribute information for that polygon.
> Does anyone have a function like that?
Bob,
The ogrinfo utility will approximately do this; however, it doesn't do
real point-in-polygon tests, it just verifies that the bounding box of
the polygon(s) intersects the target point(or region). If that is
sufficient then you would use something like:
ogrinfo -ro -spat 481710 4752312 481710 4752312 polygon.shp polygon
http://ogr.maptools.org
There is a point-in-polygon available in the CVS version of Shapelib if you
would like to construct a program to do this yourself. I haven't tried the
point in polygon code, but I assume it works ok. It may require some digging
to figure out how to use it and build an app to do the full operation.
I have attached the readme and code I received for this. It should work with
a regular release shapelib.
Best regards,
--
---------------------------------------+--------------------------------------
I set the clouds in motion - turn up | Frank Warmerdam, warmerdam@pobox.com
light and sound - activate the windows | http://pobox.com/~warmerdam
and watch the world go round - Rush | Geospatial Programmer for Rent
--------------030000040708090904020500
Content-Type: text/plain;
name="Shape_PointInPoly.cpp"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="Shape_PointInPoly.cpp"
/******************************************************************************
* $Id: Shape_PointInPoly.cpp,v 1.1 2004/01/09 16:47:57 fwarmerdam Exp $
*
* Project: Shapelib
* Purpose: Commandline program to generate points-in-polygons from a
* shapefile as a shapefile.
* Author: Marko Podgorsek, d-mon@siol.net
*
******************************************************************************
* Copyright (c) 2004, Marko Podgorsek, d-mon@siol.net
*
* This software is available under the following "MIT Style" license,
* or at the option of the licensee under the LGPL (see LICENSE.LGPL). This
* option is discussed in more detail in shapelib.html.
*
* --
*
* Permission is hereby granted, free of charge, to any person obtaining a
* copy of this software and associated documentation files (the "Software"),
* to deal in the Software without restriction, including without limitation
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
* and/or sell copies of the Software, and to permit persons to whom the
* Software is furnished to do so, subject to the following conditions:
*
* The above copyright notice and this permission notice shall be included
* in all copies or substantial portions of the Software.
*
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
* DEALINGS IN THE SOFTWARE.
******************************************************************************
*
* $Log: Shape_PointInPoly.cpp,v $
* Revision 1.1 2004/01/09 16:47:57 fwarmerdam
* New
*
*/
static char rcsid[] =
"$Id: Shape_PointInPoly.cpp,v 1.1 2004/01/09 16:47:57 fwarmerdam Exp $";
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <shapefil.h>
#define MAXINTERSECTIONPOINTS 255
enum loopDir {
kExterior,
kInterior,
kError
};
struct DPoint2d
{
DPoint2d()
{
x = y = 0.0;
};
DPoint2d(double x, double y)
{
this->x = x;
this->y = y;
};
double x,y;
};
struct IntersectPoint
{
IntersectPoint(void)
{
x = y = 0.0;
boundry_nmb = 0;
loopdirection = kError;
};
double x,y;
int boundry_nmb;
loopDir loopdirection;
};
loopDir LoopDirection(DPoint2d *vertices, int vertsize)
{
int i;
double sum = 0.0;
for(i=0;i<vertsize-1;i++)
{
sum += (vertices[i].x*vertices[i+1].y)-(vertices[i].y*vertices[i+1].x);
}
if(sum>0)
return kInterior;
else
return kExterior;
}
DPoint2d CreatePointInPoly(SHPObject *psShape, int quality)
{
int i, j, k, end, vert, pointpos;
double part, dx, xmin, xmax, ymin, ymax, y, x3, x4, y3, y4, len, maxlen = 0;
DPoint2d *vertices;
loopDir direction;
IntersectPoint mp1, mp2, point1, point2, points[MAXINTERSECTIONPOINTS];
xmin = psShape->dfXMin;
ymin = psShape->dfYMin;
xmax = psShape->dfXMax;
ymax = psShape->dfYMax;
part = (ymax-ymin)/(quality+1);
dx = xmax-xmin;
for(i=0;i<quality;i++)
{
y = ymin+part*(i+1);
pointpos = 0;
for(j=0;j<psShape->nParts;j++)
{
if(j==psShape->nParts-1)
end = psShape->nVertices;
else
end = psShape->panPartStart[j+1];
vertices = new DPoint2d [end-psShape->panPartStart[j]];
for(k=psShape->panPartStart[j],vert=0;k<end;k++)
{
vertices[vert].x = psShape->padfX[k];
vertices[vert++].y = psShape->padfY[k];
}
direction = LoopDirection(vertices, vert);
for(k=0;k<vert-1;k++)
{
y3 = vertices[k].y;
y4 = vertices[k+1].y;
if((y3 >= y && y4 < y) || (y3 <= y && y4 > y)) //I check >= only once, because if it's not checked now (y3) it will be in the next iteration (which is y4 now)
{
point1.boundry_nmb = j;
point1.loopdirection = direction;
x3 = vertices[k].x;
x4 = vertices[k+1].x;
if(y3==y)
{
point1.y = y3;
point1.x = x3;
if(direction == kInterior) //add point 2 times if the direction is interior, so that the final count of points is even
{
points[pointpos++]=point1;
}
}
else
{
point1.x = xmin+(((((x4-x3)*(y-y3))-((y4-y3)*(xmin-x3)))/((y4-y3)*dx))*dx); //striped down calculation of intersection of 2 lines
point1.y = y;
}
points[pointpos++]=point1;
}
}
delete [] vertices;
}
for(j=1;j<pointpos;j++) //sort the found intersection points by x value
{
for(k=j;k>0;k--)
{
if(points[k].x < points[k-1].x)
{
point1 = points[k];
points[k] = points[k-1];
points[k-1] = point1;
}
else
{
break;
}
}
}
for(j=0;j<pointpos-1;j++)
{
point1 = points[j];
point2 = points[j+1];
if((point1.loopdirection == kExterior && //some checkings for valid point
point2.loopdirection == kExterior &&
point1.boundry_nmb == point2.boundry_nmb &&
j%2==0) ||
(point1.loopdirection == kExterior &&
point2.loopdirection == kInterior) ||
(point1.loopdirection == kInterior &&
point2.loopdirection == kExterior))
{
len = sqrt(pow(point1.x-point2.x, 2)+pow(point1.y-point2.y, 2));
if(len >= maxlen)
{
maxlen = len;
mp1 = point1;
mp2 = point2;
}
}
}
}
return DPoint2d((mp1.x+mp2.x)*0.5, (mp1.y+mp2.y)*0.5);
}
int main(int argc, char* argv[])
{
if(argc != 3)
{
printf("Usage: %s shpfile_path quality\n", argv[0]);
return 1;
}
int i, nEntities, quality;
SHPHandle hSHP;
SHPObject *psShape;
DPoint2d pt;
quality = atoi(argv[2]);
hSHP = SHPOpen(argv[1], "rb");
SHPGetInfo(hSHP, &nEntities, NULL, NULL, NULL);
printf("PointInPoly v1.0, by Marko Podgorsek\n----------------\n");
for( i = 0; i < nEntities; i++ )
{
psShape = SHPReadObject( hSHP, i );
if(psShape->nSHPType == SHPT_POLYGON)
{
pt = CreatePointInPoly(psShape, quality);
printf("%d: x=%f y=%f\n",i, pt.x,pt.y);
}
SHPDestroyObject( psShape );
}
SHPClose(hSHP);
return 0;
}
--------------030000040708090904020500
Content-Type: text/plain;
name="Shape_PointInPoly_README.txt"
Content-Transfer-Encoding: 7bit
Content-Disposition: inline;
filename="Shape_PointInPoly_README.txt"
===============================================================================
Project: Shape_PoinInPoly
Purpose: Sample and the function for calculatin a point in a polygon
(complex,compound - it doesn't matter). Can be used for labeling.
Author: Copyright (c) 2004, Marko Podgorsek, d-mon@siol.net
===============================================================================
Requires: shapelib 1.2 (http://shapelib.maptools.org/)
Tested and created on platform:
Windows 2000 Professional
Visual Studio .NET 7.0
P4 2.4 GHz
1GB RAM
I just found out about the ShapeLib, GDAL and OGR and I must say that they're
all great projects.
I belive I'll use some of those libraries in the future. Right now I'm using
only shapelib.
The thing that led me to the http://wwww.maptools.org was the need of finding
the point in poly...but as I found out that even OGR didn't support it. So
there I was. I was forced to make my own function. Well, it was fun. I learned
a lot.
I wrote this function for the Autodesk Autocad 2004 MPolygon, because there was
no function to do this in the Object Arx SDK (the Acad programming SDK). Well,
it will be in the 2005 release...but, still. There is a function in the
Autodesk Map 2004 version...in the menu.
Not usefull when you need the coordinates, not the point on the screen...
So when the Acad version was done I was thinking of doing it on the Shape files,
too. A little bit of changing the structures and variable
types (so they're not using Object Arx) and I was done.
And here it is....Contribution from me to the ShapeLib world :)...and maybe even
OGR (a little bit of changing there).
Some statistics:
For about 69000 polygons in Autocad picture (.dwg files)
Autodesk Map 2004 was creating centroids (the menu command) about 45s (1 scan
line)
My function, with 3 scan lines took about 5s. And I was drawing the dots on the
picture...
-------------------------------------------------------------------------------
DPoint2d CreatePointInPoly(SHPObject *psShape, int quality)
The second parameter quality tell the function just how many rays shall it use
to get the point.
quality = 3 works very well, but anything below 5 is good.
This doesn't mean that the execution will slow down, but it just finds a good
point. That's all.
The qality shows on the compound objects (multiple poligons with more than one
exterior loop) - if not enough rays, then there may be no centroid.
Or the U shaped thin polygon, only the bootom (below the y center line) is fat.
Autodesk Map with one scan line will create the centroid on one of the thin
parts, because it only uses the y center line. If you have more rays, one will
surely pass the fat area and centroid will be created there.
-------------------------------------------------------------------------------
Anyone using this function:
Just send me an e-mail, so I'll see if I did anything good for the public.
And you can send me e-mail with questions also.
--------------030000040708090904020500--
This archive was generated by Pipermail. |
MapTools.org -- Hosted by DM Solutions Group |