[R-sig-Geo] 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 <vijaylulla using gmail.com>
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 <aolinto.lst using gmail.com>
> 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
>> R-sig-Geo using r-project.org
>> 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
> vlulla using iupui.edu
> ORCID: https://orcid.org/0000-0002-0823-2522
> Webpage: http://vijaylulla.com
>

	[[alternative HTML version deleted]]