8

I would like to calculate the maximum, median and minimum latidudinal coordinates for a large set of polygons, using QGIS or R. So the highest and lowest points reached by each polygon across a latidudinal range. Central coords are easy to do, but converting the polygons to points, and using a distance matrix in QGIS crashes my computer. Is there a more effecient way to do this?

Its been asked in various guises here, but either 1. Unclear 2. outdated I think version wise or 3. crashes on my PC as I have about 2000 polygons.

7

You can do it like this:

library(raster)
# example data
g <- getData('GADM', country='BRA', level=1)

ext <- t(sapply(1:length(g), function(i) as.vector(extent(g[i,]))))
colnames(ext) <- c('xmin', 'xmax', 'ymin', 'ymax')

head(ext)
#          xmin      xmax       ymin      ymax
#[1,] -73.98971 -66.58875 -11.145161 -7.121320
#[2,] -38.23634 -35.15182 -10.501529 -8.814987
#[3,] -54.87619 -49.86681  -1.236008  4.442360
#[4,] -73.79568 -56.09750  -9.814520  2.246201
#[5,] -46.61705 -37.34903 -18.349859 -8.533636
#[6,] -41.42347 -37.25208  -7.858196 -2.784583

And continue like this:

d <- data.frame(state=g$NAME_1, ext)
head(d)

#     state      xmin      xmax       ymin      ymax
#1     Acre -73.98971 -66.58875 -11.145161 -7.121320
#2  Alagoas -38.23634 -35.15182 -10.501529 -8.814987
#3    Amapá -54.87619 -49.86681  -1.236008  4.442360
#4 Amazonas -73.79568 -56.09750  -9.814520  2.246201
#5    Bahia -46.61705 -37.34903 -18.349859 -8.533636
#6    Ceará -41.42347 -37.25208  -7.858196 -2.784583
  • Nice! --------------------- – mdsumner Jan 15 '16 at 21:35
  • Well thats just elegant. All answers seem to do the trick, but this wins imo for effeciency and simplicity. Thank you! – Guest2345234562 Jan 16 '16 at 8:33
5

In QGIS, you could use the Polygon from Layer Extent... tool from the toolbar (Vector > Research Tools > Polygon from Layer Extent...).

This essentially outputs a bounding box layer for each feature (if you select the option) with fields containing the coordinates of the max, centre and min of X and Y along with few other statistics:

Polygon from Layer Extent

Polygon example

Attribute table of results

  • 1
    I think that does the trick - thank you! However, and this may be a general QGIS question - it does not conserve object ID. Is the order of the outputs the same as the input polygon? Is there a way to ensure that it is, i.e. matching each correct bounding box with the other data in the original shapefile. – Guest2345234562 Jan 14 '16 at 14:22
  • @Guest2345234562 - Most welcome buddy! Yes, I believe the order of the outputs matches the inputs. I did a quick check by spatially joining the layers (e.g. using the Join attributes by location tool) and they seem to match up nicely. – Joseph Jan 14 '16 at 14:39
  • Does anyone know where this is in QGIS 3? Thanks – Ian Allan Mar 27 at 3:54
5

In R

  1. read them with x <- gdal::readOGR(datasource, layername) from just about any format

  2. use as(x, "class") coercion to convert from polygons to their line boundaries to their component points (and handily, record the object and ring IDs)

  3. Use summary functions in the standard way for the X/Y grouped by polygon ID

Example:

library(rgdal)
dsn <- system.file("vectors/ps_cant_31.MIF", package = "rgdal")[1]
ogrInfo(dsn=dsn, layer="ps_cant_31")
ps_cant_31 <- readOGR(dsn=dsn, layer="ps_cant_31")

## cast to lines and then to points (creates columns  Lines.NR, Lines.ID,     Line.NR to identify pieces)

p <-  as(as(ps_cant_31, "SpatialLinesDataFrame"), "SpatialPointsDataFrame")
coords <- coordinates(p)
## see that Lines.NR groups your original polygons (nrow(ps_cant_31))
#spplot(p["Lines.NR"])

## summarize the Y-coordinate into groups defined by original polygon object
tapply(coords[,2], p$Lines.NR, median)
tapply(coords[,2], p$Lines.NR, max)
tapply(coords[,2], p$Lines.NR, min)

You mention "latitude range" but don't specify if your coordinates are in a projected coordinate system, if they are you need to ensure you unproject and classify by latitude (not just Y) if that is the case.

I would highly recommend you check out dplyr for the final summarizing, R's built-in tools are very powerful but tiresome in hindsight. Unfortunately with sp you need to manually convert between arrays of coordinates and data.frame form pretty routinely to get these kinds of answers, but it's all doable.

Finally, if this is actually what you are after it won't take any time at all, as long as your virtual memory resources are a match for the data set you have, to read it all in at once and do the job in one step like this.

2

In QGIS, you can use expressions to add fields using the Field Calculator

Following expressions will give you the min and max coordinates for the polygon

xmin($geometry)
xmax($geometry)
ymin($geometry)
ymax($geometry)

For center coordinates, you can use

xmin(centroid($geometry))
ymin(centroid($geometry))

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.