#Construct a global grid with cells approximately 1000 miles across
dggs <- dgconstruct(spacing=2, metric=TRUE, resround='down')
## Resolution: 15, Area (km^2): 3.55473501726709, Spacing (km): 1.86210705534755, CLS (km): 2.12744663988395
#Set working directory
#setwd("C:\\Users\\b1019175\\Desktop\\Linda\\R_skripte\\dggs-integration")
#Read the test data (point data)
sbg <- st_read("C:\\Users\\b1019175\\Desktop\\Linda\\R_skripte\\dggs-integration\\data\\alle_Wasserpunkte.shp")
## Reading layer `alle_Wasserpunkte' from data source `C:\Users\b1019175\Desktop\Linda\R_skripte\dggs-integration\data\alle_Wasserpunkte.shp' using driver `ESRI Shapefile'
## Simple feature collection with 27421 features and 14 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -1.797693e+308 ymin: -1.797693e+308 xmax: 498823 ymax: 322503.2
## Projected CRS: MGI / Austria GK M31
#st_coordinates(sbg)[,1] #show coordinates
st_is_longlat(sbg) #test whether the coordinate reference system is latlon
## [1] FALSE
st_crs(sbg) #show the coordinate reference system of the data
## Coordinate Reference System:
## User input: MGI / Austria GK M31
## wkt:
## PROJCRS["MGI / Austria GK M31",
## BASEGEOGCRS["MGI",
## DATUM["Militar-Geographische Institut",
## ELLIPSOID["Bessel 1841",6377397.155,299.1528128,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4312]],
## CONVERSION["Austria Gauss-Kruger M31",
## METHOD["Transverse Mercator",
## ID["EPSG",9807]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",13.3333333333333,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["Scale factor at natural origin",1,
## SCALEUNIT["unity",1],
## ID["EPSG",8805]],
## PARAMETER["False easting",450000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",-5000000,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["northing (X)",north,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["easting (Y)",east,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Engineering survey, topographic mapping."],
## AREA["Austria between 11°50'E and 14°50'E of Greenwich (29°30'E and 32°30'E of Ferro)."],
## BBOX[46.4,11.83,48.79,14.84]],
## ID["EPSG",31258]]
sbg_trans <- st_transform(sbg, "+proj=longlat +datum=WGS84 +no_defs") #transform coordinate reference system to latlon
st_crs(sbg_trans) #show again the crs
## Coordinate Reference System:
## User input: +proj=longlat +datum=WGS84 +no_defs
## wkt:
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]
sbg_trans <- na.omit(sbg_trans) # ignore no data values
#Make two sepparate columns for lat/lon
sbg_trans <- sbg_trans %>%
dplyr::mutate(lon = sf::st_coordinates(.)[,1],
lat = sf::st_coordinates(.)[,2])
x<-sbg_trans[1:500,]#Create a subset of 500 datapoints
#Get the corresponding grid cells for each earthquake epicenter (lat-long pair)
x$cell <- dgGEO_to_SEQNUM(dggs,x$lon,x$lat)$seqnum
#Converting SEQNUM to GEO gives the center coordinates of the cells
cellcenters <- dgSEQNUM_to_GEO(dggs,x$cell)
#Get the number of earthquakes in each cell
xcounts <- x %>% group_by(cell) %>% summarise(count=n())
#Get the grid cell boundaries for cells which had quakes
grid <- dgcellstogrid(dggs,xcounts$cell,frame=TRUE,wrapcells=TRUE)
#Update the grid cells' properties to include the number of earthquakes
#in each cell
grid <- merge(grid,xcounts,by.x="cell",by.y="cell")
#Get Salzburg admin boundaries
bl <- st_read("C:\\Users\\b1019175\\Desktop\\Linda\\R_skripte\\dggs-integration\\data\\BL_Salzburg.shp")
## Reading layer `BL_Salzburg' from data source `C:\Users\b1019175\Desktop\Linda\R_skripte\dggs-integration\data\BL_Salzburg.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 3 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 304951.8 ymin: 338346.5 xmax: 450321.9 ymax: 460283
## Projected CRS: MGI / Austria Lambert
bl <- st_transform(bl, "+proj=longlat +datum=WGS84 +no_defs")
p<-ggplot() +
geom_sf(data = bl, size = 2, shape = 23, fill = "white") +
geom_polygon(data=grid, aes(x=long, y=lat, group=group, fill = count))
p
q <-ggplot() +
geom_sf(data = bl, size = 2, shape = 23, fill = "white") +
geom_sf(data = x)
q
