[R-sig-Geo] Converting UTMs to SpatialPolygons object

[Michael Sumner]

Here we go, I'm assuming each ring is its own simple object, and a few
other things.

id = rep(seq(1,3),each=12)
az = rep(seq(0,359,30),3)

## assume you meant to use all ids here?
d = data.frame(id,az)
d$dist = round(rnorm(36,10,3),1)
xp = function(azimuth, r){cos(((90-azimuth)*pi)/180)*r}
yp =  function(azimuth, r){sin(((90-azimuth)*pi)/180)*r}
ctr_e = rnorm(3,667314,50)
ctr_n = rnorm(3,4784480,50)

d$utm_e = round(xp(d$az,d$dist),1) + ctr_e[d$id]
d$utm_n = round(yp(d$az,d$dist),1) + ctr_n[d$id]

library(sp)
myCRS = CRS("+proj=utm +zone=18 +ellps=WGS84")

## create list of coords, including repeat of first row for each
l.coords <- lapply(split(d, d$id), function(x)
as.matrix(x[c(seq_len(nrow(x)), 1), c("utm_e", "utm_n")]))


l.Polygons <- vector("list",length(l.coords))

for (i in seq_len(length(l.coords))) {
    l.Polygons[[i]] <- Polygons(list(Polygon(l.coords[i])), as.character(i))
}

SpPolys <- SpatialPolygons(l.Polygons, proj4string = myCRS)

plot(SpPolys, col = c("grey", "red", "blue"))


On Mon, Jul 16, 2012 at 10:24 AM,  <seth at swbigelow.net> wrote:
>
>
>         It seems my data-generating code was cut off. Here's another try:
>
>         id = rep(seq(1,3),each=12)            # Generate polygon
> IDs'
> az = rep(seq(0,359,30),3)               # Generate
> azimuth sequence
> c = data.frame(id
> [1],az)                        # Combine in
> data frame
> c$dist = round(rnorm(36,10,3),1)  # Generate random distances from
> center of polygon
>
>         xp = function(azimuth, r){cos(((90-azimuth)*pi)/180)*r}
>       # Function to calculate horizontal distance from polygon
> center
>
>         yp =  function(azimuth, r){sin(((90-azimuth)*pi)/180)*r}
>       # Function to calculate vertical distance from polygon
> center
>
>         ctr_e = rnorm(3,667314,50)       # Randomly select eastings
> for polygon centers
> ctr_n = rnorm(3,4784480,50)     # Randomly select northings for
> polygon centers
>
>         c$utm_e = round(xp(c$az,c$dist),1) + ctr_e[c$id]  # Calculate
> point eastings
> c$utm_n = round(yp(c$az,c$dist),1) + ctr_n[c$id]  # Calculate point
> northings
>
>         myCRS = CRS("+proj=utm +zone=18 +ellps=WGS84")    # define
> proj4string
>
>         Appreciatively,
>
>         Seth Bigelow
>
>
>
> Links:
> ------
> [1] http://sitemail.hostway.com/http:/
>
>
>         [[alternative HTML version deleted]]
>
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>



-- 
Michael Sumner
Hobart, Australia
e-mail: mdsumner at gmail.com