I am trying to extend GADM (admin level 2) districts out into the sea using Thiessen polygons (or voronoi I think if you are a mathematician) up to a boundary e.g. the EEZ zones/200km, and then be able to merge these to create 'coastal zones' incorporating both the admin district and the corresponding coastal waters. I've found a way to do this in ArcPy online, but I'm trying to develop using open source software. I've come up with the following based on different forums and responses, and it seems to work, but the polygons seem to not line up with the original admin districts that well (maybe because I had to simplify them due to lack of processing power) and the shapefile just looks a bit odd. If anyone can immediately flag up any issues or suggest any better ways to go about this that would be super helpful! Thank you!
#Create regions; script extends admin areas to include associated coastal waters up to a boundary (200km or EEZ)
library(raster)
library(sf)
#import shapefiles of admin district and EEZ
library(rgdal)
library(rmapshaper)
library(maptools)
GADM_2<-readOGR(dsn="/Users/amycampbell/Documents/YGT_Project/Data/Regions/", layer="GADM_2")
GADM = ms_simplify(GADM_2, keep = 0.01) #this is the proportion of points to maintain- ideally want higher but then takes much longer (don;t have high processing power available)
zone<- readOGR(dsn="/Users/amycampbell/Documents/YGT_Project/Data/Regions/", layer="India_200kmline") #also have an EEZ shapefile which could be used instead, but 200km buffer was used to keep it simpler for now
#code adapted from https://stackoverflow.com/questions/37738251/how-to-expand-the-polygon-to-reach-nearby-line-in-r
#create nearest neighbour polygons
#get polygon coordinates
library(dismo)
g<-unique(geom(GADM))
print(g)
#Make thiessen polygons up to the edge of the zone
v<-voronoi(g[,c('x', 'y')], ext=extent(zone))
plot(v)
#assign names from GADM to new polygons
v$name<-GADM$NAME_2[v$id]
#aggregate (dissolve) polygons by the name
a<- aggregate(v, 'name')
#assign projection to these polygons
proj4string(a)<-CRS("+proj=lcc +lat_1=12.472955 +lat_2=35.17280444444444 +lat_0=24 +lon_0=80 +x_0=4000000 +y_0=4000000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")
a2<-spTransform(v, CRS("+proj=lcc +lat_1=12.472955 +lat_2=35.17280444444444 +lat_0=24 +lon_0=80 +x_0=4000000 +y_0=4000000 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0"))
#remove areas outside of zone
i<-crop(a2,zone)
plot(i)
#save as shapefile (both methods seem to work!)
writeOGR(i, dsn="/Users/amycampbell/Documents/YGT_Project/Data/Regions/", layer="india_regions", driver="ESRI Shapefile")
#shapefile(i, filename="/Users/amycampbell/Documents/YGT_Project/Data/Regions/india_regions.shp", overwrite=TRUE)
)! 1 nm is 1.852 km.