How Can I Generate a Triangular Network with R?
Part of My Coding Diary with R
In this article, we will demonstrate how you can generate a triangular network with the famous open source software R.
Why Do I Need a Triangular Network?
As I have said in my introduction article, I am planning to build a functional principal component analysis on a set of data.
While that does not mean that I will necessary need a triangular network, there are some unique properties I would like to use later. The most important properties of a functional EOF set based on triangular network is that it is continuous and the cross product is easy to compute.
Furthermore, if I want to develop a basic finite element method in the future, I probably will have to generate a triangular network anyway.
The Packages I Used in This Article
library(“tripack”, lib.loc=”~/R/win-library/3.4")
library(“maptools”, lib.loc=”~/R/win-library/3.4")
library(“raster”, lib.loc=”~/R/win-library/3.4")
library(“rgeos”, lib.loc=”~/R/win-library/3.4")The package “tripack” was used for triangular network generation.
The package “maptools” was used to read the shape file of Taiwan.
“raster” and “rgeos” were used to perform most of the geo-statistical operations such as distance counting or combining polygons. I can’t remember which functions belong to which package specifically.
Anyways, let’s get started.
Reading a Shape File in R
So I downloaded a shape file of Taiwan and found the function to read it.
Reading a shape file in R wasn’t difficult at all. And combining all the polygons into one greater raster object wasn’t hard either. The functions readShapePoly() and aggregate() did their job well.
taiwan=readShapePoly(“COUNTY_MOI_1060525.shp”)
taiwan.union <- aggregate(taiwan)The real problem was to know how R read these shape file. It was stored into a raster object. It had a very confusing data structure:
It took me some time to figure out the right way to call the coordinates of a specific polygon.
The coordinate of the k-th point in the j-th polygon of the i-th polygon set in the raster object obj would be:
obj@polygons[[i]]@Polygons[[j]]@coords[k,]After applying the raster object to aggregate(), there will only be one polygon set in that object.
However, since there are a lot of islands, I had to find the polygon with the most points to know which of the polygon was the main island.
counter=c()
max=1
prev=1
for(i in 1:length(taiwan.union@polygons[[1]]@Polygons)){
max=max(c(length(taiwan.union@polygons[[1]]@Polygons[[i]]@coords[,1]),max))
if(max>prev){
counter=i
}
prev=max
}
plot(taiwan.union@polygons[[1]]@Polygons[[counter]]@coords,asp=1,type=”l”)The plot() function on the last line above verified that this searching rationale is valid.
I actually didn’t need a resolution that high, so I reduced the points on the polygon to an amount of 500:
L=length(taiwan.union@polygons[[1]]@Polygons[[counter]]@coords)/2
border=matrix(nrow=500,ncol=2)
border=taiwan.union@polygons[[1]]@Polygons[[counter]]@coords[L%/%500*seq(1,500),]Generating Points within the Polygon Using Geo-distance Contours
To generate a triangular network, you have to have some points within the polygon. I generated the points using contour lines of the geo-distance to the shore.
Unfortunately the gDistance() function only returns the value when the input point is outside the input polygon, so it took me some time to find an algorithm to generate the complement of Taiwan island.
xlim=c(120,122)
ylim=c(21.5,25.5)
x_in=seq(xlim[1],xlim[2],by=.1)
y_in=seq(ylim[1],ylim[2],by=.1)land <- raster(extent(list(x=xlim,y=ylim)))
res(land) <- .1
proj4string(land)<-proj4string(taiwan.union)
landpolygon <- rasterToPolygons(land)
land_com <- gDifference(landpolygon,taiwan.union)
Basically the code above is just telling the computer to build a rectangular raster object, and then subtract it with the raster object I had before.
In the graph below, you can see this is a complement to the original polygon by the small gap on the northeastern edge.
Moving on, I generated a regular grid and calculate the distance of each grid point to the shore:
dis=matrix(nrow=length(x_in),ncol=length(y_in))
for(i in 1:length(x_in)){
for(j in 1:length(y_in)){
dis[i,j]=gDistance(readWKT(paste(“POINT”,”(“,xlim[1]+(i-1)*.1,” “,ylim[1]+(j-1)*.1,”)”,sep=””)),land_com)
}
}And then I generated the corresponding geo-distance contour lines:
dis.contour=contourLines(dis,levels=seq(0,.6,.05))
for(i in 1:length(dis.contour)){
dis.contour[[i]]$x=xlim[1]+dis.contour[[i]]$x*(xlim[2]-xlim[1])
dis.contour[[i]]$y=ylim[1]+dis.contour[[i]]$y*(ylim[2]-ylim[1])
}Plotting everything together confirmed the validity of our coding so far.
plot(border,asp=1,type=”l”)
for(i in 1:length(dis.contour)){
lines(dis.contour[[i]]$x,dis.contour[[i]]$y,text(dis.contour[[i]]$x[5*i],dis.contour[[i]]$y[5*i],labels=dis.contour[[i]]$level),type=”l”)
}
Generating the Triangular Network and Delete Irrelevant Triangles
To generate the triangular network, I first had to combine all the points together into one matrix:
points=border
for(i in 1:length(dis.contour)){
dis.matrix=cbind(dis.contour[[i]]$x[1:(length(dis.contour[[i]]$x)-1)],dis.contour[[i]]$y[1:(length(dis.contour[[i]]$y)-1)])
points=rbind(points,dis.matrix)
}Then I called the tri.mesh() function to generate a triangular network with the matrix. This process was very straightforward.
tin_Taiwan=tri.mesh(x=points[,1],y=points[,2])
plot.tri(tin_Taiwan,asp=1)
lines(taiwan.union,col="red")But when I plotted the triangular network, something very wrong immediately appeared on the screen.
I had to come up with an algorithm to delete the triangles that were outside the border of the original polygon.
So I turned to another function, triangles(), which turns a triangulation object to a matrix containing the same information.
The matrix would be much easier to compute with, because each row represents a triangle, and the first three columns represents the three points of a specific triangle.
Therefore, I could filter out the triangles with all three border points. This algorithm would also delete some of the triangle on the edges that were still inside the island, but the overall error would still improve.
tin_Taiwan_matrix=triangles(tin_Taiwan)
for(i in 1:nrow(tin_Taiwan_matrix)){
if(tin_Taiwan_matrix[i,1]<500&&tin_Taiwan_matrix[i,2]<500&&tin_Taiwan_matrix[i,3]<500){
tin_Taiwan_matrix[i,]=NA
}
}And that was basically it, though to plot the resulting network also required some effort.
plot(c(points[tin_Taiwan_matrix[1,1],1],points[tin_Taiwan_matrix[1,2],1]),c(points[tin_Taiwan_matrix[1,1],2],points[tin_Taiwan_matrix[1,2],2]),type=”l”,xlim=c(120,122),ylim=c(21.5,25.5),asp=1)
lines(c(points[tin_Taiwan_matrix[1,2],1],points[tin_Taiwan_matrix[1,3],1]),c(points[tin_Taiwan_matrix[1,2],2],points[tin_Taiwan_matrix[1,3],2]))
lines(c(points[tin_Taiwan_matrix[1,3],1],points[tin_Taiwan_matrix[1,1],1]),c(points[tin_Taiwan_matrix[1,3],2],points[tin_Taiwan_matrix[1,1],2]))for(i in 2:nrow(tin_Taiwan_matrix)){
if(!is.na(tin_Taiwan_matrix[i,1])){
lines(c(points[tin_Taiwan_matrix[i,1],1],points[tin_Taiwan_matrix[i,2],1]),c(points[tin_Taiwan_matrix[i,1],2],points[tin_Taiwan_matrix[i,2],2]))
lines(c(points[tin_Taiwan_matrix[i,2],1],points[tin_Taiwan_matrix[i,3],1]),c(points[tin_Taiwan_matrix[i,2],2],points[tin_Taiwan_matrix[i,3],2]))
lines(c(points[tin_Taiwan_matrix[i,3],1],points[tin_Taiwan_matrix[i,1],1]),c(points[tin_Taiwan_matrix[i,3],2],points[tin_Taiwan_matrix[i,1],2]))
}
}
lines(taiwan.union,col=”red”)
And at last I got the graph which I have shown at the very beginning of this article:
What is Next?
Now that I get my triangular network, I would like to generate a functional principal component analysis with this network. I will explain what this really means and how I deal with it.
About This Series…
This article is part of the series My Coding Diary with R. To know how the series has been developing, check the intro article out:
