#Include libraries
library(dggridR)
library(collapse)

#Construct a global grid with cells approximately 1000 miles across
dggs          <- dgconstruct(spacing=1000, metric=FALSE, resround='down')

#Load included test data set
data(dgquakes)

#Get the corresponding grid cells for each earthquake epicenter (lat-long pair)
dgquakes$cell <- dgGEO_to_SEQNUM(dggs,dgquakes$lon,dgquakes$lat)$seqnum

#Converting SEQNUM to GEO gives the center coordinates of the cells
cellcenters   <- dgSEQNUM_to_GEO(dggs,dgquakes$cell)

#Get the number of earthquakes in each cell
quakecounts   <- dgquakes |> fcount(cell, name = "count")

#Get the grid cell boundaries for cells which had quakes
grid          <- dgcellstogrid(dggs,quakecounts$cell)

#Update the grid cells' properties to include the number of earthquakes
#in each cell
grid          <- merge(grid,quakecounts,by.x="seqnum",by.y="cell")

#Make adjustments so the output is more visually interesting
grid$count    <- log(grid$count)
cutoff        <- fquantile(grid$count, 0.9)
grid          <- grid |> fmutate(count = ifelse(count>cutoff,cutoff,count))

#Get polygons for each country of the world
countries <- map_data("world")

#Plot everything on a flat map

# Handle cells that cross 180 degrees
wrapped_grid = st_wrap_dateline(grid, options = c("WRAPDATELINE=YES","DATELINEOFFSET=180"), quiet = TRUE)

ggplot() +
    geom_polygon(data=countries, aes(x=long, y=lat, group=group), fill=NA, color="black")   +
    geom_sf     (data=wrapped_grid, aes(fill=count), color=alpha("white", 0.4)) +
    geom_point  (aes(x=cellcenters$lon_deg, y=cellcenters$lat_deg)) +
    scale_fill_gradient(low="blue", high="red")

library(sf)

#Get the grid cell boundaries for the whole Earth using this dggs in a form
#suitable for printing to a KML file
grid <- dgearthgrid(dggs)

#Update the grid cells' properties to include the number of earthquakes
#in each cell
grid$count <- merge(grid, quakecounts, by.x="seqnum", by.y="cell", all.x=TRUE)

#Write out the grid
st_write(grid, "quakes_per_cell.kml", layer="quakes", driver="KML")

#Include libraries
library(dggridR)

N <- 100    #How many cells to sample

#Distribute the points uniformly on a sphere using equations from
#http://mathworld.wolfram.com/SpherePointPicking.html
u     <- runif(N)
v     <- runif(N)
theta <- 2*pi*u      * 180/pi
phi   <- acos(2*v-1) * 180/pi
lon   <- theta-180
lat   <- phi-90

df    <- data.frame(lat=lat,lon=lon)

#Construct a global grid in which every hexagonal cell has an area of
#100,000 miles^2. You could, of course, choose a much smaller value, but these
#will show up when I map them later.

#Note: Cells can only have certain areas, the `dgconstruct()` function below
#will tell you which area is closest to the one you want. You can also round
#up or down.

#Note: 12 cells are actually pentagons with an area 5/6 that of the hexagons
#But, with millions and millions of hexes, you are unlikely to choose one
#Future versions of the package will make it easier to reject the pentagons
dggs    <- dgconstruct(area=100000, metric=FALSE, resround='nearest')

#Get the corresponding grid cells for each randomly chosen lat-long
df$cell <- dgGEO_to_SEQNUM(dggs,df$lon,df$lat)$seqnum

#Get the hexes for each of these cells
gridfilename <- dgcellstogrid(dggs,df$cell)

#Get the grid in a more convenient format
grid <- dgcellstogrid(dggs,df$cell)
grid <- st_wrap_dateline(grid, options = c("WRAPDATELINE=YES","DATELINEOFFSET=180"), quiet = TRUE)

#Get polygons for each country of the world
countries <- map_data("world")

#Plot everything on a flat map
p <- ggplot() +
    geom_polygon(data=countries, aes(x=long, y=lat, group=group), fill=NA, color="black")   +
    geom_sf(data=grid, fill=alpha("green", alpha=0.4), color=alpha("white", alpha=0.4))
p

#Include libraries
library(dggridR)

N <- 100    #How many cells to sample

#Construct a global grid in which every hexagonal cell has an area of
#100,000 miles^2. You could, of course, choose a much smaller value, but these
#will show up when I map them later.

#Note: Cells can only have certain areas, the `dgconstruct()` function below
#will tell you which area is closest to the one you want. You can also round
#up or down.

#Note: 12 cells are actually pentagons with an area 5/6 that of the hexagons
#But, with millions and millions of hexes, you are unlikely to choose one
#Future versions of the package will make it easier to reject the pentagons
dggs    <- dgconstruct(area=100000, metric=FALSE, resround='nearest')

maxcell <- dgmaxcell(dggs)                     #Get maximum cell id
cells   <- sample(1:maxcell, N, replace=FALSE) #Choose random cells
grid    <- dgcellstogrid(dggs,cells)           #Get grid

#Get the grid in a more convenient format
grid <- dgcellstogrid(dggs,df$cell)
grid <- st_wrap_dateline(grid, options = c("WRAPDATELINE=YES","DATELINEOFFSET=180"), quiet = TRUE)

#Get polygons for each country of the world
countries <- map_data("world")

#Plot everything on a flat map
p <- ggplot() +
    geom_polygon(data=countries, aes(x=long, y=lat, group=group), fill=NA, color="black")   +
    geom_sf(data=grid, fill=alpha("green", 0.4), color=alpha("white", 0.4))
p

library(dggridR)
#Generate a global grid whose cells are ~100,000 miles^2
dggs         <- dgconstruct(area=100000, metric=FALSE, resround='nearest')
#Save the cells to a KML file for use in other software
gridfilename <- dgearthgrid(dggs,savegrid=tempfile())

library(dggridR)

#Generate a dggs specifying an intercell spacing of ~25 miles
dggs      <- dgconstruct(spacing=100, metric=FALSE, resround='nearest')

#Read in the South Africa's borders from the shapefile
sa_border <- st_read(dg_shpfname_south_africa(), layer="ZAF_adm0")
st_crs(sa_border) = 4326

#Get a grid covering South Africa
sa_grid   <- dgshptogrid(dggs, dg_shpfname_south_africa())

#Plot South Africa's borders and the associated grid
p <- ggplot() +
    geom_sf(data=sa_border, fill=NA, color="black")   +
    geom_sf(data=sa_grid, fill=alpha("blue", 0.4), color=alpha("white", 0.4))
p