# [R-sig-Geo] From spatial points to "spatial pictures"

boB Rudis bob at rudis.net
Wed Aug 5 17:33:28 CEST 2015

```Maurizio provided his data set offline. The code + data + sample plots
are in this gist:
https://gist.github.com/hrbrmstr/e29cf9138e480db9e67d

The data is just as he described it, so it's really just a matter of
creating some polygons from the the points spline interpolation:

library(ggplot2)
library(ggthemes)
library(sp)

#' Make smooth polys from points
#'
#' From: http://gis.stackexchange.com/a/24929/29544
#'
#' @param xy n x 2 matrix with n >= k
#' @param k # of vertices to expand points to
#' @param ... other args passed to \code{spline()}
spline.poly <- function(xy, vertices, k=3, ...) {

# Wrap k vertices around each end.
n <- dim(xy)
if (k >= 1) {
data <- rbind(xy[(n-k+1):n,], xy, xy[1:k, ])
} else {
data <- xy
}

# Spline the x and y coordinates.
data.spline <- spline(1:(n+2*k), data[,1], n=vertices, ...)
x <- data.spline\$x
x1 <- data.spline\$y
x2 <- spline(1:(n+2*k), data[,2], n=vertices, ...)\$y

# Retain only the middle part.
cbind(x1, x2)[k < x & x <= n+k, ]
}

# read in the data
trees <- read.csv("CrownsExample.csv")

# sort data by tree & point
trees <- trees[with(trees, order(TREE, POINT)),]

# plot just points colored by tree to make sure points look ok
ggplot(trees, aes(x=X, y=Y, color=factor(TREE))) + coord_equal() +
geom_point()

# small function to avoid anonymous fnction in `by` make a polygon
from our points
make_tree_poly <- function(tree) {
tree_pts <- spline.poly(as.matrix(tree[, c("X", "Y")]), 100)
Polygons(list(Polygon(tree_pts)), ID=unique(tree\$TREE))
}

# for each tree, make a set of smoothed polygons then make the whole thing
# a SpatialPolygonsDataFrame so more spatial ops can be done on it (or more
# data stored with it)
tree_polys <- SpatialPolygonsDataFrame(
SpatialPolygons(unclass(by(trees, trees\$TREE, make_tree_poly)),
proj4string=CRS("+proj=utm +zone=32 +ellps=GRS80
+units=m +no_defs ")),
data.frame(id=unique(trees\$TREE), row.names=unique(trees\$TREE)))

# plot the new polygon tree crowns colored by tree
tree_map <- fortify(tree_polys)
ggplot() +
geom_map(data=tree_map, map=tree_map,
aes(x=long, y=lat, fill=id, map_id=id, color=id)) +
geom_point(data=trees, aes(x=X, y=Y)) + # original points just
to show it's accurate
coord_equal() +
theme_map()

On Wed, Aug 5, 2015 at 9:56 AM, Maurizio Marchi
<mauriziomarchi85 at gmail.com> wrote:
> Hi everybody,
> I'm maurizio Marchi and I'm working as a post-doc at the Forestry Research
> Centre of Arezzo (Italy) and I have a question concerning spatial points
> management.
> During my research activities I have to collect data about the crown
> projection on the ground of trees. The aim is to "measure" the crown area
> and to draw a sort of picture on GIS software (QGIS in my case) to analyze
> the area and the perimeter of the crowns compared to the "free area" of a
> plot (generally a circular plot) to have an idea of how much photosynthetic
> radiation is able to reach the soil. My working steps are:
> 1) During the survey on the field each crown of each tree is initially
> "drawn" collecting 8 points that have to be connected to become a spatial
> polygon:
>
> ###---example of table to be processed
> TREE   POINT   X   Y
>     1            1       x    y
>     1            2       x    y
>     1            3       x    y
>     1            4       x    y
>     1            5       x    y
>     1            6       x    y
>     1            7       x    y
>     1            8       x    y
>     2            1       x    y
>     2            2       x    y
>     2            3       x    y
>     2            4       x    y
>     2            5       x    y
>     2            6       x    y
>     2            7       x    y
>     2            8       x    y
> [......................................]
>     N            8       x    y
>
> 2) the points,sorted by Tree become the vertices of the polygons
>
> 3) The polygons have to be processed to obtain a smoothed polygons,
> circumscribing the polygonal crowns (e.g.
> http://grasswiki.osgeo.org/wiki/V.generalize_tutorial)
>
> I generally do it on QGIS and GRASS but is there something more appropriate
> on R to do this? for example with rgdal or sp package?
>
> Many thanks
>
>
> --
> Maurizio Marchi
> Calenzano (FI) - Italy
> ID Skype: maurizioxyz
> Ubuntu 14.04 LTS
> linux user 552742
>
>         [[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

```

More information about the R-sig-Geo mailing list