12

I have a village polygon layer which has more than 6,00,000 records. I have calculated centroid of each village. I want to find the distance between centroid and the farthest node of each polygon. Check the image below for reference. Black lines are polygon boundaries. enter image description here

enter image description here

enter image description here

  • interesting... I just did this friday with postgis to produce a circle around a polygon. I need a few minutes to look for the code I used.. i.stack.imgur.com/EKnkg.png – kttii Jun 27 '16 at 13:57
  • 1
    First we might need to know wich programs you have at your disposal. Also how did you create thoses centroids and nodes? (Even if it seems a bit obvious that the nodes on the polygons are the ones used to set the boundaries of your shapes, but did you add aditionnal point on top of thoses?) – Moreau Colin Jun 27 '16 at 13:58
  • Is the location of the centroid important? How did you create them? – GISGe Jun 27 '16 at 14:08
  • Possible duplicate - gis.stackexchange.com/questions/133099/… – klewis Jun 27 '16 at 14:27
  • If the centroid truly is central, then it's the radius of the smallest circle centred on that point that fits the polygon (en.wikipedia.org/wiki/Smallest-circle_problem) – Mark Ireland Jun 27 '16 at 15:04
3

It looks like you are working with MapInfo. The "Distance Calculator" (find it in tool manager) has many options and I think it can handle this task too. There is an article about it here: http://web.pb.com/mapinfopro-jul-2013/Toolbox-Distance-Calculator

| improve this answer | |
15

Using PostGIS, I used ST_ConvexHull to simplify the polygon for a faster result:

Get the furthest Point:

SELECT Villages_v4_Trial_region.geom as FarPoint from (
SELECT ST_PointN(ST_ExteriorRing(ST_ConvexHull(Villages_v4_Trial_region.geom)),
generate_series(1, ST_NPoints(ST_ExteriorRing(ST_ConvexHull(Villages_v4_Trial_region.geom))))) as points, 
geom
FROM Villages_v4_Trial_region
ORDER BY ST_MaxDistance(points,ST_Centroid(Villages_v4_Trial_region.geom)) DESC
LIMIT 1;

And if you are interested in creating a Circle from the centroid:

SELECT ST_Buffer(Center,ST_Distance(Center,FarPoint)) as Circle
FROM (
SELECT Villages_v4_Trial_region.geom as FarPoint, Center from (
    SELECT ST_PointN(ST_ExteriorRing(ST_ConvexHull(Villages_v4_Trial_region.geom)),
    generate_series(1, ST_NPoints(ST_ExteriorRing(ST_ConvexHull(Villages_v4_Trial_region.geom))))) as points,
    ST_Centroid(Villages_v4_Trial_region.geom) as Center, 
    geom
    FROM Villages_v4_Trial_region
    ) as Villages_v4_Trial_region
    ORDER BY ST_MaxDistance(points,Center) DESC
    LIMIT 1) as foo;

enter image description here

| improve this answer | |
  • Simple, fast, efficient. Thank you for posting this because this'll help me too in what i am working on right now. – Moreau Colin Jun 27 '16 at 14:39
  • @kttii I don't know how to use PostGIS. Can you provide a simpler solution, in arc or mapinfo or qgis – Divya Jun 28 '16 at 4:33
  • @kttii So I installed Postgresql. I copy-pasted this exact query but it gave an ERROR: column "the_geom" does not exist. What do I do? – Divya Jun 28 '16 at 9:34
  • the_geom should be replaced by your geometry field name. You would have to put your data into PostgreSQL as well. PostgreSQL is a database like MSSQL. PostGIS is an extension to make the database spatially aware and provide all the ST_ functions. – kttii Jun 28 '16 at 12:52
  • @kttii I updated the field name from the_geom to "gid" in my database. After running the query again, I got this ERROR: function st_convexhull(integer) does not exist – Divya Jun 29 '16 at 7:25
4

Using next PyQGIS code:

from math import sqrt

layer = iface.activeLayer()

feats = [ feat for feat in layer.getFeatures() ]

n = len(feats)

centroids = [ feat.geometry().centroid().asPoint() for feat in feats ]
polygons = [ feat.geometry().asPolygon()[0] for feat in feats ]

lengths = []

for i, pol in enumerate(polygons):
    max_dist = 0
    idx_j = 0
    for j, point in enumerate(pol):
        dist = sqrt(centroids[i].sqrDist(point))
        if dist > max_dist:
            max_dist = dist
            idx_j = j
    print i, idx_j, max_dist
    lengths.append([centroids[i], pol[idx_j]])

crs = layer.crs()
epsg = crs.postgisSrid()

uri = "LineString?crs=epsg:" + str(epsg) + "&field=id:integer""&index=yes"

mem_layer = QgsVectorLayer(uri,
                           'max_distance',
                           'memory')

prov = mem_layer.dataProvider()

feats = [ QgsFeature() for i in range(n) ]

for i, feat in enumerate(feats):
    feat.setAttributes([i])
    feat.setGeometry(QgsGeometry.fromPolyline(lengths[i]))

prov.addFeatures(feats)

QgsMapLayerRegistry.instance().addMapLayer(mem_layer)

and this shapefile (with 11 features):

enter image description here

I got a memory layer where polylines were the distance between centroid and farthest point of each polygon (feature); as it can be observed at the next image:

enter image description here

At the Python Console of QGIS, it was also printed the index of feature, the index of point in feature where distance from centroid is a maximum and, finally, maximum distance.

enter image description here

| improve this answer | |
  • I don't know how to use PyQGIS. Can you provide a simpler solution, in arc or mapinfo or qgis – Divya Jun 28 '16 at 4:34
  • 1
    Try this link for help with getting started with PyQgis spatialgalaxy.net/2014/10/09/… – kttii Jun 28 '16 at 13:02
0

As it looks like you're using MapInfo, here is a MapBasic function I wrote a while back for an in-house tool I was working on. It takes a source node (your centroid point) and a region object (the polygon) as arguments and returns a point object at the furthest node in the polygon from the source point.

Function GetFurthest(ByVal oNode1 as Object, ByVal oObj as Object) as Object

Dim sourceE,sourceN,East,North,Longest,Dist as Float,
    nNodes,nPolys,i,j as SmallInt,
    oNode2 as Object

    sourceE = CentroidX(oNode1)
    sourceN = CentroidY(oNode1)
    Longest = 0

    nPolys = ObjectInfo(oObj,OBJ_INFO_NPOLYGONS)
    For i = 1 to nPolys
        nNodes = ObjectInfo(oObj,OBJ_INFO_NPOLYGONS+nPolys)
        For j = 1 to nNodes
            East = ObjectNodeX(oObj,i,j)
            North = ObjectNodeY(oObj,i,j)
            Dist = Distance(sourceE,sourceN,East,North,"m")
            If Dist > Longest then
                Longest = Dist
                oNode2 = CreatePoint(East,North)
            End if
        Next
    Next

    GetFurthest = oNode2

End Function
| improve this answer | |

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.