read them with
x <- gdal::readOGR(datasource, layername) from just about any format
as(x, "class") coercion to convert from polygons to their line boundaries to their component points (and handily, record the object and ring IDs)
Use summary functions in the standard way for the X/Y grouped by polygon ID
dsn <- system.file("vectors/ps_cant_31.MIF", package = "rgdal")
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))
## 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.