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 |
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 |
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 |
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 |
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 |
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 |
Free forum by Nabble | Edit this page |