How Can I Generate a Triangular Network with R?

Why Do I Need a Triangular Network?

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")

Reading a Shape File in R

taiwan=readShapePoly(“COUNTY_MOI_1060525.shp”)
taiwan.union <- aggregate(taiwan)
obj@polygons[[i]]@Polygons[[j]]@coords[k,]
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”)
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

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)
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)
}
}
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])
}
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

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)
}
tin_Taiwan=tri.mesh(x=points[,1],y=points[,2])
plot.tri(tin_Taiwan,asp=1)
lines(taiwan.union,col="red")
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
}
}
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”)

What is Next?

About This Series…

A Taiwanese student who studied Renewable Energy in Freiburg. Now studying smart distribution grids / energy systems in Trondheim.