drawing a polygon from several lines

classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

drawing a polygon from several lines

Antonio Silva
Hello all,

I want to select the cells of a grid that in a range of latitudes and
depths.
Latitudes limits are indicated by two lines and the depths by two isobaths
from a shapefile.
The shapefile can be downloaded at
https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0

Below are the steps I've already taken. I could not find a way to create a
polygon from the lines of depth and latitude. With this polygon I could
easily select the overlapped cells.

I would really appreciate any help. Best regards

----
library(rgeos)
library(sp)
rm(list = ls())

# import shape
isobaths <- readOGR(".","isobath")
isobaths <- spTransform(bath, CRS("+proj=longlat +datum=WGS84"))
proj4string(isobaths)
summary(isobaths)

# create a spatial grid and polygons
grd <-
GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
polys <- as.SpatialPolygons.GridTopology(grd)
proj4string(polys) <- CRS("+proj=longlat +datum=WGS84")
summary(polys)
IDs <- sapply(slot(polys, "polygons"), function(x) slot(x, "ID"))

# create spatial lines
# southern limit
x <- c(-47.5,-45.5)
y <- c(-24.8,-24.8)
ls <- SpatialLines(list(Lines(Line(cbind(x,y)),
ID="b")),proj4string=CRS("+proj=longlat +datum=WGS84"))

# northern limit
x <- c(-46.2,-44.5)
y <- c(-24.02,-24.02)
ln <- SpatialLines(list(Lines(Line(cbind(x,y)),
ID="c")),proj4string=CRS("+proj=longlat +datum=WGS84"))

# plot isolines polygons and lines
plot(polys)
axis(1);axis(4)
text(coordinates(polys), labels=IDs, cex=0.7)
plot(isobaths,add=T)
plot(ls,add=T,col="blue",lwd=2)
plot(ln,add=T,col="orange",lwd=2)

# select isobaths
summary(isobaths)
dmin <- subset(isobaths, ID %in% "-25")
proj4string(dmin) <- CRS("+proj=longlat +datum=WGS84")

dmax <- subset(isobaths, ID %in% "-60")
proj4string(dmax) <- CRS("+proj=longlat +datum=WGS84")

plot(dmin,add=T,col="red",lwd=2)
plot(dmax,add=T,col="red",lwd=2)

# from here I'd like to have a polygon to select the cells.
na.omit(over(polys,dmin))

--
Antonio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

        [[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: drawing a polygon from several lines

Antonio Silva
Dear colleagues

I finally could select the squares (cells) under the polygon. I extract the
points from the selected isobaths that fell within lat / lon limits and
with them I created a spatial polygon.

Bellow I share the script I wrote. Considering that I don't know much about
mapping in R, I would appreciate to hear suggestions to improve the code.

Best regards.

--
Antônio Olinto Ávila da Silva
Instituto de Pesca (Fisheries Institute)
São Paulo, Brasil

===================

library(rgdal)
library(rgeos)
rm(list = ls())

# import shapes
# isobath
https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0
isobaths <- readOGR(".","isobath")
isobaths <- spTransform(isobaths, CRS("+proj=longlat +datum=WGS84"))
proj4string(isobaths)
summary(isobaths)

# create a spatial grid and polygons
grd <-
GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
squares <- as.SpatialPolygons.GridTopology(grd)
proj4string(squares) <- CRS("+proj=longlat +datum=WGS84")
summary(squares)
IDs <- sapply(slot(squares, "polygons"), function(x) slot(x, "ID"))

# set limits
latmin <- -24.8
latmax <- -24.02
lonmin <- -47.2
lonmax <- -44.8
depmin <- -25
depmax <- -60

# bounding box
mat.area <-
matrix(c(lonmin,latmax,lonmax,latmax,lonmax,latmin,lonmin,latmin),ncol=2,byrow=T)
spp.area <- SpatialPolygons(list(Polygons(list(Polygon(mat.area)) ,
ID='1')))
proj4string(spp.area) <- CRS("+proj=longlat +datum=WGS84")

# plot polygons isolines and bounding box
plot(squares, axes=T)
#~ text(coordinates(squares),labels=IDs, cex=0.7)
plot(isobaths,add=T)
plot(spp.area,border="red",add=T)

# select isobaths
summary(isobaths)
isomin <- subset(isobaths, ID %in% depmin)
proj4string(isomin) <- CRS("+proj=longlat +datum=WGS84")
isomax <- subset(isobaths, ID %in% depmax)
proj4string(isomax) <- CRS("+proj=longlat +datum=WGS84")

plot(isomin,add=T,col="deepskyblue3",lwd=2)
plot(isomax,add=T,col="dodgerblue4",lwd=2)

# cut isobaths
limmin <- subset(coordinates(isomin)[[1]][[1]],
          coordinates(isomin)[[1]][[1]][,2]<=latmax &
coordinates(isomin)[[1]][[1]][,2]>=latmin &
          coordinates(isomin)[[1]][[1]][,1]<=lonmax &
coordinates(isomin)[[1]][[1]][,1]>=lonmin)
limmax <- subset(coordinates(isomax)[[1]][[1]],
          coordinates(isomax)[[1]][[1]][,2]<=latmax &
coordinates(isomax)[[1]][[1]][,2]>=latmin &
          coordinates(isomax)[[1]][[1]][,1]<=lonmax &
coordinates(isomax)[[1]][[1]][,1]>=lonmin)

points(limmin,col="chartreuse3")
points(limmax,col="chartreuse4")

# create the polygon
spp.farea <-
SpatialPolygons(list(Polygons(list(Polygon(rbind(limmin,limmax))),ID='1')))
proj4string(spp.farea) <- CRS("+proj=longlat +datum=WGS84")

plot(spp.farea,col="chocolate3",lwd=2,add=T)

# select the squares under the polygon
fareasq <- squares[spp.farea,]
fareace <- gCentroid(fareasq,byid=TRUE)

# final plot
plot(fareasq,axes=T)
points(fareace,pch=10,col="darkgreen",cex=4)
text(coordinates(squares), labels=IDs, cex=0.7)
plot(isomin,add=T,col="deepskyblue3",lwd=2)
plot(isomax,add=T,col="dodgerblue4",lwd=2)
plot(spp.area,border="red",add=T,lty=2)

        [[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: drawing a polygon from several lines

Vijay Lulla
Very nice, Antonio!!  Thanks for sharing your script.  The only minor
suggestion I have to offer is that you can make your "cut isobath" portion
easier to read (at least IMO) by using `with` and `between` (either
dplyr::between or data.table::between).  Below is how I would've done it.

limmin <- with(list(coords = coordinates(isomin)[[1]][[1]]),
               subset(coords, between(coords[, 1], lonmin, lonmax) &
                              between(coords[, 2], latmin, latmax)))

Again, thanks for sharing your problem and solution.  I found it
instructive!
Cordially,
Vijay.


On Fri, Aug 24, 2018 at 9:16 AM Antonio Silva <[hidden email]> wrote:

> Dear colleagues
>
> I finally could select the squares (cells) under the polygon. I extract the
> points from the selected isobaths that fell within lat / lon limits and
> with them I created a spatial polygon.
>
> Bellow I share the script I wrote. Considering that I don't know much about
> mapping in R, I would appreciate to hear suggestions to improve the code.
>
> Best regards.
>
> --
> Antônio Olinto Ávila da Silva
> Instituto de Pesca (Fisheries Institute)
> São Paulo, Brasil
>
> ===================
>
> library(rgdal)
> library(rgeos)
> rm(list = ls())
>
> # import shapes
> # isobath
> https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0
> isobaths <- readOGR(".","isobath")
> isobaths <- spTransform(isobaths, CRS("+proj=longlat +datum=WGS84"))
> proj4string(isobaths)
> summary(isobaths)
>
> # create a spatial grid and polygons
> grd <-
>
> GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
> squares <- as.SpatialPolygons.GridTopology(grd)
> proj4string(squares) <- CRS("+proj=longlat +datum=WGS84")
> summary(squares)
> IDs <- sapply(slot(squares, "polygons"), function(x) slot(x, "ID"))
>
> # set limits
> latmin <- -24.8
> latmax <- -24.02
> lonmin <- -47.2
> lonmax <- -44.8
> depmin <- -25
> depmax <- -60
>
> # bounding box
> mat.area <-
>
> matrix(c(lonmin,latmax,lonmax,latmax,lonmax,latmin,lonmin,latmin),ncol=2,byrow=T)
> spp.area <- SpatialPolygons(list(Polygons(list(Polygon(mat.area)) ,
> ID='1')))
> proj4string(spp.area) <- CRS("+proj=longlat +datum=WGS84")
>
> # plot polygons isolines and bounding box
> plot(squares, axes=T)
> #~ text(coordinates(squares),labels=IDs, cex=0.7)
> plot(isobaths,add=T)
> plot(spp.area,border="red",add=T)
>
> # select isobaths
> summary(isobaths)
> isomin <- subset(isobaths, ID %in% depmin)
> proj4string(isomin) <- CRS("+proj=longlat +datum=WGS84")
> isomax <- subset(isobaths, ID %in% depmax)
> proj4string(isomax) <- CRS("+proj=longlat +datum=WGS84")
>
> plot(isomin,add=T,col="deepskyblue3",lwd=2)
> plot(isomax,add=T,col="dodgerblue4",lwd=2)
>
> # cut isobaths
> limmin <- subset(coordinates(isomin)[[1]][[1]],
>           coordinates(isomin)[[1]][[1]][,2]<=latmax &
> coordinates(isomin)[[1]][[1]][,2]>=latmin &
>           coordinates(isomin)[[1]][[1]][,1]<=lonmax &
> coordinates(isomin)[[1]][[1]][,1]>=lonmin)
> limmax <- subset(coordinates(isomax)[[1]][[1]],
>           coordinates(isomax)[[1]][[1]][,2]<=latmax &
> coordinates(isomax)[[1]][[1]][,2]>=latmin &
>           coordinates(isomax)[[1]][[1]][,1]<=lonmax &
> coordinates(isomax)[[1]][[1]][,1]>=lonmin)
>
> points(limmin,col="chartreuse3")
> points(limmax,col="chartreuse4")
>
> # create the polygon
> spp.farea <-
> SpatialPolygons(list(Polygons(list(Polygon(rbind(limmin,limmax))),ID='1')))
> proj4string(spp.farea) <- CRS("+proj=longlat +datum=WGS84")
>
> plot(spp.farea,col="chocolate3",lwd=2,add=T)
>
> # select the squares under the polygon
> fareasq <- squares[spp.farea,]
> fareace <- gCentroid(fareasq,byid=TRUE)
>
> # final plot
> plot(fareasq,axes=T)
> points(fareace,pch=10,col="darkgreen",cex=4)
> text(coordinates(squares), labels=IDs, cex=0.7)
> plot(isomin,add=T,col="deepskyblue3",lwd=2)
> plot(isomax,add=T,col="dodgerblue4",lwd=2)
> plot(spp.area,border="red",add=T,lty=2)
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> [hidden email]
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>


--
Vijay Lulla

Assistant Professor,
Dept. of Geography, IUPUI
425 University Blvd, CA-207C.
Indianapolis, IN-46202
[hidden email]
ORCID: https://orcid.org/0000-0002-0823-2522
Webpage: http://vijaylulla.com

        [[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: drawing a polygon from several lines

Antonio Silva
It's my pleasure to share Lulla. Thanks for the suggestion, I will adopt it.

Regards,

Antonio

Em sex, 24 de ago de 2018 às 10:58, Vijay Lulla <[hidden email]>
escreveu:

> Very nice, Antonio!!  Thanks for sharing your script.  The only minor
> suggestion I have to offer is that you can make your "cut isobath" portion
> easier to read (at least IMO) by using `with` and `between` (either
> dplyr::between or data.table::between).  Below is how I would've done it.
>
> limmin <- with(list(coords = coordinates(isomin)[[1]][[1]]),
>                subset(coords, between(coords[, 1], lonmin, lonmax) &
>                               between(coords[, 2], latmin, latmax)))
>
> Again, thanks for sharing your problem and solution.  I found it
> instructive!
> Cordially,
> Vijay.
>
>
> On Fri, Aug 24, 2018 at 9:16 AM Antonio Silva <[hidden email]>
> wrote:
>
>> Dear colleagues
>>
>> I finally could select the squares (cells) under the polygon. I extract
>> the
>> points from the selected isobaths that fell within lat / lon limits and
>> with them I created a spatial polygon.
>>
>> Bellow I share the script I wrote. Considering that I don't know much
>> about
>> mapping in R, I would appreciate to hear suggestions to improve the code.
>>
>> Best regards.
>>
>> --
>> Antônio Olinto Ávila da Silva
>> Instituto de Pesca (Fisheries Institute)
>> São Paulo, Brasil
>>
>> ===================
>>
>> library(rgdal)
>> library(rgeos)
>> rm(list = ls())
>>
>> # import shapes
>> # isobath
>> https://www.dropbox.com/sh/us859uyijns36r9/AAC9N4OSwweLuW-fdsBWUy-ra?dl=0
>> isobaths <- readOGR(".","isobath")
>> isobaths <- spTransform(isobaths, CRS("+proj=longlat +datum=WGS84"))
>> proj4string(isobaths)
>> summary(isobaths)
>>
>> # create a spatial grid and polygons
>> grd <-
>>
>> GridTopology(cellcentre.offset=c(-47.75,-25.416667),cellsize=c(10/60,10/60),cells.dim=c(23,12))
>> squares <- as.SpatialPolygons.GridTopology(grd)
>> proj4string(squares) <- CRS("+proj=longlat +datum=WGS84")
>> summary(squares)
>> IDs <- sapply(slot(squares, "polygons"), function(x) slot(x, "ID"))
>>
>> # set limits
>> latmin <- -24.8
>> latmax <- -24.02
>> lonmin <- -47.2
>> lonmax <- -44.8
>> depmin <- -25
>> depmax <- -60
>>
>> # bounding box
>> mat.area <-
>>
>> matrix(c(lonmin,latmax,lonmax,latmax,lonmax,latmin,lonmin,latmin),ncol=2,byrow=T)
>> spp.area <- SpatialPolygons(list(Polygons(list(Polygon(mat.area)) ,
>> ID='1')))
>> proj4string(spp.area) <- CRS("+proj=longlat +datum=WGS84")
>>
>> # plot polygons isolines and bounding box
>> plot(squares, axes=T)
>> #~ text(coordinates(squares),labels=IDs, cex=0.7)
>> plot(isobaths,add=T)
>> plot(spp.area,border="red",add=T)
>>
>> # select isobaths
>> summary(isobaths)
>> isomin <- subset(isobaths, ID %in% depmin)
>> proj4string(isomin) <- CRS("+proj=longlat +datum=WGS84")
>> isomax <- subset(isobaths, ID %in% depmax)
>> proj4string(isomax) <- CRS("+proj=longlat +datum=WGS84")
>>
>> plot(isomin,add=T,col="deepskyblue3",lwd=2)
>> plot(isomax,add=T,col="dodgerblue4",lwd=2)
>>
>> # cut isobaths
>> limmin <- subset(coordinates(isomin)[[1]][[1]],
>>           coordinates(isomin)[[1]][[1]][,2]<=latmax &
>> coordinates(isomin)[[1]][[1]][,2]>=latmin &
>>           coordinates(isomin)[[1]][[1]][,1]<=lonmax &
>> coordinates(isomin)[[1]][[1]][,1]>=lonmin)
>> limmax <- subset(coordinates(isomax)[[1]][[1]],
>>           coordinates(isomax)[[1]][[1]][,2]<=latmax &
>> coordinates(isomax)[[1]][[1]][,2]>=latmin &
>>           coordinates(isomax)[[1]][[1]][,1]<=lonmax &
>> coordinates(isomax)[[1]][[1]][,1]>=lonmin)
>>
>> points(limmin,col="chartreuse3")
>> points(limmax,col="chartreuse4")
>>
>> # create the polygon
>> spp.farea <-
>>
>> SpatialPolygons(list(Polygons(list(Polygon(rbind(limmin,limmax))),ID='1')))
>> proj4string(spp.farea) <- CRS("+proj=longlat +datum=WGS84")
>>
>> plot(spp.farea,col="chocolate3",lwd=2,add=T)
>>
>> # select the squares under the polygon
>> fareasq <- squares[spp.farea,]
>> fareace <- gCentroid(fareasq,byid=TRUE)
>>
>> # final plot
>> plot(fareasq,axes=T)
>> points(fareace,pch=10,col="darkgreen",cex=4)
>> text(coordinates(squares), labels=IDs, cex=0.7)
>> plot(isomin,add=T,col="deepskyblue3",lwd=2)
>> plot(isomax,add=T,col="dodgerblue4",lwd=2)
>> plot(spp.area,border="red",add=T,lty=2)
>>
>>         [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> [hidden email]
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>>
>
>
> --
> Vijay Lulla
>
> Assistant Professor,
> Dept. of Geography, IUPUI
> 425 University Blvd, CA-207C.
> Indianapolis, IN-46202
> [hidden email]
> ORCID: https://orcid.org/0000-0002-0823-2522
> Webpage: http://vijaylulla.com
>

        [[alternative HTML version deleted]]

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