# Why Do I Need a Triangular Network?

`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=1prev=1for(i in 1:length(taiwan.union@polygons[]@Polygons)){ max=max(c(length(taiwan.union@polygons[]@Polygons[[i]]@coords[,1]),max)) if(max>prev){ counter=i } prev=max}plot(taiwan.union@polygons[]@Polygons[[counter]]@coords,asp=1,type=”l”)`
`L=length(taiwan.union@polygons[]@Polygons[[counter]]@coords)/2border=matrix(nrow=500,ncol=2)border=taiwan.union@polygons[]@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,xlim,by=.1)y_in=seq(ylim,ylim,by=.1)land <- raster(extent(list(x=xlim,y=ylim)))res(land) <- .1proj4string(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+(i-1)*.1,” “,ylim+(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+dis.contour[[i]]\$x*(xlim-xlim) dis.contour[[i]]\$y=ylim+dis.contour[[i]]\$y*(ylim-ylim)}`
`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=borderfor(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”)`