Polygon width

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

Polygon width

Paulo Flores Ribeiro
Hello,

I have a shapefile with ca. 25000 polygons. Each polygon has an average
of 40 vertices (nodes). I would like to extract, for each polygon, the
distance separating the two most distant vertices (aka "polygon
diagonal" or "maximum polygon width"). It is not important whether the
polygon is convex or concave, so the lines connecting the vertices can
be inside or outside the polygon. The desired result would be a
two-column array, with a number of rows equal to the number of polygons,
and where the first column is the id of the polygons, and the second the
"maximum width" of the corresponding polygon.

What would be the best way to do this, considering that the calculation
will probably require frequent updates (e.g. due to changes in the shape
of the polygons)?

Thanks in advance,

PauloFR

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Polygon width

Barry Rowlingson
Do you want great-circle distance or is your space small enough that you
can use planar coordinates?

Are your polygons all single rings or are there islands and/or holes? Does
that matter?

The straightforward way would be to coerce the polygons to points, compute
the distance matrix, then take the maximum. Depending on if you are reading
your shapefile into sp or sf classes, the functions would be a bit
different. You should try and implement the straightforward way, test it
for correctness, and then worry about the "best way" if the straightforward
way isn't what you need. Its often the case that "best" ways need fancy
data structures or complex algorithms with more opportunity for bugs. Start
simple, work up.

For example, using sf, here's the max distance between any points in the
first feature of `pcs`

> max(st_distance(st_cast(st_geometry(pcs[1,]),"POINT")))
172.556 m

loop from 1 to N or otherwise apply over the features, and you're done.

Barry



On Wed, Apr 25, 2018 at 11:26 AM, Paulo Flores Ribeiro <
[hidden email]> wrote:

> Hello,
>
> I have a shapefile with ca. 25000 polygons. Each polygon has an average of
> 40 vertices (nodes). I would like to extract, for each polygon, the
> distance separating the two most distant vertices (aka "polygon diagonal"
> or "maximum polygon width"). It is not important whether the polygon is
> convex or concave, so the lines connecting the vertices can be inside or
> outside the polygon. The desired result would be a two-column array, with a
> number of rows equal to the number of polygons, and where the first column
> is the id of the polygons, and the second the "maximum width" of the
> corresponding polygon.
>
> What would be the best way to do this, considering that the calculation
> will probably require frequent updates (e.g. due to changes in the shape of
> the polygons)?
>
> Thanks in advance,
>
> PauloFR
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Polygon width

Paulo Flores Ribeiro
Thank you, Barry. I am using planar coordinates (meters) and Euclidean
distances. The shape or size of the polygons is not important (for
illustration purposes, imagine that polygons are farm boundaries and
that each farm is a single polygon). I know how to extract the maximum
distance in a single polygon by calculating the distance matrix and
selecting the maximum value. My "difficulty" is exactly in the coding of
the loop process (or in the construction of the "function" to be used in
the apply approach), so that I can apply it "automatically" to 25,000
polygons. I am using the rgadl package, but I can switch to sf.
Thanks for any help.
Cheers,
PauloFR

Às 12:27 de 25-04-2018, Barry Rowlingson escreveu:

> Do you want great-circle distance or is your space small enough that
> you can use planar coordinates?
>
> Are your polygons all single rings or are there islands and/or holes?
> Does that matter?
>
> The straightforward way would be to coerce the polygons to points,
> compute the distance matrix, then take the maximum. Depending on if
> you are reading your shapefile into sp or sf classes, the functions
> would be a bit different. You should try and implement the
> straightforward way, test it for correctness, and then worry about the
> "best way" if the straightforward way isn't what you need. Its often
> the case that "best" ways need fancy data structures or complex
> algorithms with more opportunity for bugs. Start simple, work up.
>
> For example, using sf, here's the max distance between any points in
> the first feature of `pcs`
>
> > max(st_distance(st_cast(st_geometry(pcs[1,]),"POINT")))
> 172.556 m
>
> loop from 1 to N or otherwise apply over the features, and you're done.
>
> Barry
>
>
>
> On Wed, Apr 25, 2018 at 11:26 AM, Paulo Flores Ribeiro
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     Hello,
>
>     I have a shapefile with ca. 25000 polygons. Each polygon has an
>     average of 40 vertices (nodes). I would like to extract, for each
>     polygon, the distance separating the two most distant vertices
>     (aka "polygon diagonal" or "maximum polygon width"). It is not
>     important whether the polygon is convex or concave, so the lines
>     connecting the vertices can be inside or outside the polygon. The
>     desired result would be a two-column array, with a number of rows
>     equal to the number of polygons, and where the first column is the
>     id of the polygons, and the second the "maximum width" of the
>     corresponding polygon.
>
>     What would be the best way to do this, considering that the
>     calculation will probably require frequent updates (e.g. due to
>     changes in the shape of the polygons)?
>
>     Thanks in advance,
>
>     PauloFR
>
>     _______________________________________________
>     R-sig-Geo mailing list
>     [hidden email] <mailto:[hidden email]>
>     https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>     <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>
>
>
> <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
> Sem vírus. www.avg.com
> <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
>
>
> <#DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>


        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Polygon width

Barry Rowlingson
In reply to this post by Barry Rowlingson
Loop over the row indexes of an sf-class object `pctest`, applying a
function given the row index which returns the distance:

 > ppa =
sapply(1:nrow(pctest),function(i){max(st_distance(st_cast(st_geometry(pctest[i,]),"POINT")))})
 > ppa
  [1] 172.55598 360.77081 107.53889 137.17785  51.66645 132.82052 113.00875
  [8] 161.02432 141.13909  88.67002

If you want to make that a new column, do `pctest$maxdist = ppa`

And test on a simple example - make sure a unit square returns
approximately sqrt(2)!

Barry




On Wed, Apr 25, 2018 at 5:24 PM, Paulo Flores Ribeiro <
[hidden email]> wrote:

> Thank you, Barry. I am using planar coordinates (meters) and Euclidean
> distances. The shape or size of the polygons is not important (for
> illustration purposes, imagine that polygons are farm boundaries and that
> each farm is a single polygon). I know how to extract the maximum
> distance in a single polygon by calculating the distance matrix and
> selecting the maximum value. My "difficulty" is exactly in the coding of
> the loop process (or in the construction of the "function" to be used in
> the apply approach), so that I can apply it "automatically" to 25,000
> polygons. I am using the rgadl package, but I can switch to sf.
> Thanks for any help.
> Cheers,
> PauloFR
>
> Às 12:27 de 25-04-2018, Barry Rowlingson escreveu:
>
> Do you want great-circle distance or is your space small enough that you
> can use planar coordinates?
>
> Are your polygons all single rings or are there islands and/or holes? Does
> that matter?
>
> The straightforward way would be to coerce the polygons to points, compute
> the distance matrix, then take the maximum. Depending on if you are reading
> your shapefile into sp or sf classes, the functions would be a bit
> different. You should try and implement the straightforward way, test it
> for correctness, and then worry about the "best way" if the straightforward
> way isn't what you need. Its often the case that "best" ways need fancy
> data structures or complex algorithms with more opportunity for bugs. Start
> simple, work up.
>
> For example, using sf, here's the max distance between any points in the
> first feature of `pcs`
>
> > max(st_distance(st_cast(st_geometry(pcs[1,]),"POINT")))
> 172.556 m
>
> loop from 1 to N or otherwise apply over the features, and you're done.
>
> Barry
>
>
>
> On Wed, Apr 25, 2018 at 11:26 AM, Paulo Flores Ribeiro <
> [hidden email]> wrote:
>
>> Hello,
>>
>> I have a shapefile with ca. 25000 polygons. Each polygon has an average
>> of 40 vertices (nodes). I would like to extract, for each polygon, the
>> distance separating the two most distant vertices (aka "polygon diagonal"
>> or "maximum polygon width"). It is not important whether the polygon is
>> convex or concave, so the lines connecting the vertices can be inside or
>> outside the polygon. The desired result would be a two-column array, with a
>> number of rows equal to the number of polygons, and where the first column
>> is the id of the polygons, and the second the "maximum width" of the
>> corresponding polygon.
>>
>> What would be the best way to do this, considering that the calculation
>> will probably require frequent updates (e.g. due to changes in the shape of
>> the polygons)?
>>
>> Thanks in advance,
>>
>> PauloFR
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
>
> <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient> Sem
> vírus. www.avg.com
> <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
> <#m_-2085622191602335703_DAB4FAD8-2DD7-40BB-A1B8-4E2AA1F9FDF2>
>
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Polygon width

Paulo Flores Ribeiro
Thank you, Barry! In the meanwhile, I've been exploring other paths and
I think I found another solution, using loop. For a
SpatialPolygonsDataFrame named "map4" the code was:

dmax <- data.frame()
for (i in 1:nrow(map4)) {
   poly = map4@polygons[[i]]@Polygons[[1]]@coords
   dmax <- rbind.data.frame(dmax, max(dist(cbind(poly[,1],poly[,2]))))
}
print(dmax)

Anyway, I'll certainly try out your alternative solution, based on the
(s)apply function, to see which one is most efficient.

Thank you so much again,

PauloFR


Às 18:16 de 25-04-2018, Barry Rowlingson escreveu:

> Loop over the row indexes of an sf-class object `pctest`, applying a
> function given the row index which returns the distance:
>
>  > ppa =
> sapply(1:nrow(pctest),function(i){max(st_distance(st_cast(st_geometry(pctest[i,]),"POINT")))})
>  > ppa
>   [1] 172.55598 360.77081 107.53889 137.17785  51.66645 132.82052
> 113.00875
>   [8] 161.02432 141.13909  88.67002
>
> If you want to make that a new column, do `pctest$maxdist = ppa`
>
> And test on a simple example - make sure a unit square returns
> approximately sqrt(2)!
>
> Barry
>
>
>
>
> On Wed, Apr 25, 2018 at 5:24 PM, Paulo Flores Ribeiro
> <[hidden email] <mailto:[hidden email]>> wrote:
>
>     Thank you, Barry. I am using planar coordinates (meters) and
>     Euclidean distances. The shape or size of the polygons is not
>     important (for illustration purposes, imagine that polygons are
>     farm boundaries and that each farm is a single polygon). I know
>     how to extract the maximum distance in a single polygon by
>     calculating the distance matrix and selecting the maximum value.
>     My "difficulty" is exactly in the coding of the loop process (or
>     in the construction of the "function" to be used in the apply
>     approach), so that I can apply it "automatically" to 25,000
>     polygons. I am using the rgadl package, but I can switch to sf.
>     Thanks for any help.
>     Cheers,
>     PauloFR
>
>     Às 12:27 de 25-04-2018, Barry Rowlingson escreveu:
>>     Do you want great-circle distance or is your space small enough
>>     that you can use planar coordinates?
>>
>>     Are your polygons all single rings or are there islands and/or
>>     holes? Does that matter?
>>
>>     The straightforward way would be to coerce the polygons to
>>     points, compute the distance matrix, then take the maximum.
>>     Depending on if you are reading your shapefile into sp or sf
>>     classes, the functions would be a bit different. You should try
>>     and implement the straightforward way, test it for correctness,
>>     and then worry about the "best way" if the straightforward way
>>     isn't what you need. Its often the case that "best" ways need
>>     fancy data structures or complex algorithms with more opportunity
>>     for bugs. Start simple, work up.
>>
>>     For example, using sf, here's the max distance between any points
>>     in the first feature of `pcs`
>>
>>     > max(st_distance(st_cast(st_geometry(pcs[1,]),"POINT")))
>>     172.556 m
>>
>>     loop from 1 to N or otherwise apply over the features, and you're
>>     done.
>>
>>     Barry
>>
>>
>>
>>     On Wed, Apr 25, 2018 at 11:26 AM, Paulo Flores Ribeiro
>>     <[hidden email]
>>     <mailto:[hidden email]>> wrote:
>>
>>         Hello,
>>
>>         I have a shapefile with ca. 25000 polygons. Each polygon has
>>         an average of 40 vertices (nodes). I would like to extract,
>>         for each polygon, the distance separating the two most
>>         distant vertices (aka "polygon diagonal" or "maximum polygon
>>         width"). It is not important whether the polygon is convex or
>>         concave, so the lines connecting the vertices can be inside
>>         or outside the polygon. The desired result would be a
>>         two-column array, with a number of rows equal to the number
>>         of polygons, and where the first column is the id of the
>>         polygons, and the second the "maximum width" of the
>>         corresponding polygon.
>>
>>         What would be the best way to do this, considering that the
>>         calculation will probably require frequent updates (e.g. due
>>         to changes in the shape of the polygons)?
>>
>>         Thanks in advance,
>>
>>         PauloFR
>>
>>         _______________________________________________
>>         R-sig-Geo mailing list
>>         [hidden email] <mailto:[hidden email]>
>>         https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>         <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
>>
>>
>>
>>     <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
>>     Sem vírus. www.avg.com
>>     <http://www.avg.com/email-signature?utm_medium=email&utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
>>
>>
>
>

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
Reply | Threaded
Open this post in threaded view
|

Re: Polygon width

Barry Rowlingson
In reply to this post by Barry Rowlingson
Growing a data frame is usually a bad idea, and R can get slow as the
system may have to create a fresh data frame each time. If you want to use
a for-loop, make a vector or data frame of the right size beforehand and
then fill it element-by-element - sketch:

n = nrow(data)
result = rep(NA, n)
for(i in 1:n){
  result[i] = something(i)
}




On Wed, Apr 25, 2018 at 6:31 PM, Paulo Flores Ribeiro <
[hidden email]> wrote:

> Thank you, Barry! In the meanwhile, I've been exploring other paths and
> I think I found another solution, using loop. For a
> SpatialPolygonsDataFrame named "map4" the code was:
>
> dmax <- data.frame()
> for (i in 1:nrow(map4)) {
>    poly = map4@polygons[[i]]@Polygons[[1]]@coords
>    dmax <- rbind.data.frame(dmax, max(dist(cbind(poly[,1],poly[,2]))))
> }
> print(dmax)
>
> Anyway, I'll certainly try out your alternative solution, based on the
> (s)apply function, to see which one is most efficient.
>
> Thank you so much again,
>
> PauloFR
>
>
> Às 18:16 de 25-04-2018, Barry Rowlingson escreveu:
> > Loop over the row indexes of an sf-class object `pctest`, applying a
> > function given the row index which returns the distance:
> >
> >  > ppa =
> > sapply(1:nrow(pctest),function(i){max(st_distance(
> st_cast(st_geometry(pctest[i,]),"POINT")))})
> >  > ppa
> >   [1] 172.55598 360.77081 107.53889 137.17785  51.66645 132.82052
> > 113.00875
> >   [8] 161.02432 141.13909  88.67002
> >
> > If you want to make that a new column, do `pctest$maxdist = ppa`
> >
> > And test on a simple example - make sure a unit square returns
> > approximately sqrt(2)!
> >
> > Barry
> >
> >
> >
> >
> > On Wed, Apr 25, 2018 at 5:24 PM, Paulo Flores Ribeiro
> > <[hidden email] <mailto:[hidden email]>>
> wrote:
> >
> >     Thank you, Barry. I am using planar coordinates (meters) and
> >     Euclidean distances. The shape or size of the polygons is not
> >     important (for illustration purposes, imagine that polygons are
> >     farm boundaries and that each farm is a single polygon). I know
> >     how to extract the maximum distance in a single polygon by
> >     calculating the distance matrix and selecting the maximum value.
> >     My "difficulty" is exactly in the coding of the loop process (or
> >     in the construction of the "function" to be used in the apply
> >     approach), so that I can apply it "automatically" to 25,000
> >     polygons. I am using the rgadl package, but I can switch to sf.
> >     Thanks for any help.
> >     Cheers,
> >     PauloFR
> >
> >     Às 12:27 de 25-04-2018, Barry Rowlingson escreveu:
> >>     Do you want great-circle distance or is your space small enough
> >>     that you can use planar coordinates?
> >>
> >>     Are your polygons all single rings or are there islands and/or
> >>     holes? Does that matter?
> >>
> >>     The straightforward way would be to coerce the polygons to
> >>     points, compute the distance matrix, then take the maximum.
> >>     Depending on if you are reading your shapefile into sp or sf
> >>     classes, the functions would be a bit different. You should try
> >>     and implement the straightforward way, test it for correctness,
> >>     and then worry about the "best way" if the straightforward way
> >>     isn't what you need. Its often the case that "best" ways need
> >>     fancy data structures or complex algorithms with more opportunity
> >>     for bugs. Start simple, work up.
> >>
> >>     For example, using sf, here's the max distance between any points
> >>     in the first feature of `pcs`
> >>
> >>     > max(st_distance(st_cast(st_geometry(pcs[1,]),"POINT")))
> >>     172.556 m
> >>
> >>     loop from 1 to N or otherwise apply over the features, and you're
> >>     done.
> >>
> >>     Barry
> >>
> >>
> >>
> >>     On Wed, Apr 25, 2018 at 11:26 AM, Paulo Flores Ribeiro
> >>     <[hidden email]
> >>     <mailto:[hidden email]>> wrote:
> >>
> >>         Hello,
> >>
> >>         I have a shapefile with ca. 25000 polygons. Each polygon has
> >>         an average of 40 vertices (nodes). I would like to extract,
> >>         for each polygon, the distance separating the two most
> >>         distant vertices (aka "polygon diagonal" or "maximum polygon
> >>         width"). It is not important whether the polygon is convex or
> >>         concave, so the lines connecting the vertices can be inside
> >>         or outside the polygon. The desired result would be a
> >>         two-column array, with a number of rows equal to the number
> >>         of polygons, and where the first column is the id of the
> >>         polygons, and the second the "maximum width" of the
> >>         corresponding polygon.
> >>
> >>         What would be the best way to do this, considering that the
> >>         calculation will probably require frequent updates (e.g. due
> >>         to changes in the shape of the polygons)?
> >>
> >>         Thanks in advance,
> >>
> >>         PauloFR
> >>
> >>         _______________________________________________
> >>         R-sig-Geo mailing list
> >>         [hidden email] <mailto:[hidden email]>
> >>         https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>         <https://stat.ethz.ch/mailman/listinfo/r-sig-geo>
> >>
> >>
> >>
> >>     <http://www.avg.com/email-signature?utm_medium=email&
> utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
> >>      Sem vírus. www.avg.com
> >>     <http://www.avg.com/email-signature?utm_medium=email&
> utm_source=link&utm_campaign=sig-email&utm_content=emailclient>
> >>
> >>
> >
> >
>
>

        [[alternative HTML version deleted]]

_______________________________________________
R-sig-Geo mailing list
[hidden email]
https://stat.ethz.ch/mailman/listinfo/r-sig-geo