Create mesh in R for kriging in gstat
lat long
7.16 124.21
8.6 123.35
8.43 124.28
8.15 125.08
Consider these coordinates, these coordinates correspond to weather stations that measure rain data.
The gstat input in R uses the meuse dataset. At some point in this tutorial: https://rpubs.com/nabilabd/118172 , guys are using "meuse.grid" in this line of code:
data("meuse.grid")
I don't have such a file and I don't know how to create it, can I create it using these coordinates? Or at least point me to a material that discusses how to create a custom grid for a custom area (i.e. not use administrative boundaries from GADM).
Maybe this is wrong, I don't even know if this question makes sense to R savvy people. However, I would like to hear some direction, or at least some advice. Many thanks!
Total noob in R and statistics.
EDIT: See the sample table I wrote in the tutorial I want to do.
EDIT 2: Would this method be viable? https://rstudio-pubs-static.s3.amazonaws.com/46259_d328295794034414944deea60552a942.html
source to share
I'm going to share my approach to create a kriging mesh. There are probably more efficient or graceful ways to accomplish the same task, but I hope this is a starting point for discussion.
The original poster was thinking about 1 km for every 10 pixels, but that's probably too much. I'm going to create a mesh with a mesh size of 1 km * 1 km. Also, the original poster does not specify the origin of the grid, so I'll take some time to figure out a good starting point. I also assume that the Spherical Mercator projection coordinate system is a suitable choice for projection. This is a general forecast for Google Maps or Open Street Maps.
1. Download packages
I am going to use the following packages. sp
, rgdal
and raster
are packages that provide many useful spatial analysis functions. leaflet
and mapview
are packages for fast exploratory spatial data visualization.
# Load packages
library(sp)
library(rgdal)
library(raster)
library(leaflet)
library(mapview)
2. Exploratory visualization of station locations
I created an interactive map to check the location of the four stations. Since the original poster provided the latitude and longitude of these four stations, I can create SpatialPointsDataFrame
with a Lat / Long projection . Please note that EPSG code is for latitude / longitude forecast 4326
. To learn more about EPSG code, see this tutorial ( https://www.nceas.ucsb.edu/~frazier/RSpatialGuides/OverviewCoordinateReferenceSystems.pdf ).
# Create a data frame showing the **Latitude/Longitude**
station <- data.frame(lat = c(7.16, 8.6, 8.43, 8.15),
long = c(124.21, 123.35, 124.28, 125.08),
station = 1:4)
# Convert to SpatialPointsDataFrame
coordinates(station) <- ~long + lat
# Set the projection. They were latitude and longitude, so use WGS84 long-lat projection
proj4string(station) <- CRS("+init=epsg:4326")
# View the station location using the mapview function
mapview(station)
The function mapview
will create an interactive map. We can use this map to determine what might be appropriate for the mesh to occur.
3. Determine the beginning
After checking the map, I figured the origin might be near longitude 123
and latitude 7
. This origin will be in the lower left corner of the grid. Now I need to find a coordinate representing the same point in the projection of the spherical Mercator.
# Set the origin
ori <- SpatialPoints(cbind(123, 7), proj4string = CRS("+init=epsg:4326"))
# Convert the projection of ori
# Use EPSG: 3857 (Spherical Mercator)
ori_t <- spTransform(ori, CRSobj = CRS("+init=epsg:3857"))
I first created an object SpatialPoints
based on the latitude and longitude of the origin. After that I used spTransform
to perform the conversion of the project. The object is ori_t
now a spherical Mercator projection source. Please note that the EPSG code for Spherical Mercator 3857
.
To see the value of the coordinates, we can use the function coordinates
as follows.
coordinates(ori_t)
coords.x1 coords.x2
[1,] 13692297 781182.2
4. Determine the degree of mesh
Now I need to define the extent of the mesh that can cover all four points and the desired area for kriging, which depends on the cell size and the number of cells. The following code sets the scope based on the information. I figured the cell size is 1 km * 1 km, but I need to experiment with what would be a good cell number for both the x and y direction.
# The origin has been rounded to the nearest 100
x_ori <- round(coordinates(ori_t)[1, 1]/100) * 100
y_ori <- round(coordinates(ori_t)[1, 2]/100) * 100
# Define how many cells for x and y axis
x_cell <- 250
y_cell <- 200
# Define the resolution to be 1000 meters
cell_size <- 1000
# Create the extent
ext <- extent(x_ori, x_ori + (x_cell * cell_size), y_ori, y_ori + (y_cell * cell_size))
Depending on how much I created, I can create a raster layer with a number equal to 0
. Then I can use the function again mapview
to see if the raster and four stations fit well.
# Initialize a raster layer
ras <- raster(ext)
# Set the resolution to be
res(ras) <- c(cell_size, cell_size)
ras[] <- 0
# Project the raster
projection(ras) <- CRS("+init=epsg:3857")
# Create interactive map
mapview(station) + mapview(ras)
I repeated this process several times. Finally, I decided that the number of cells is 250
both 200
for the x- and y-directions, respectively.
5. Creating a spatial mesh
Now I have created a bitmap layer with the proper degree. I can first save this bitmap as GeoTiff for future use.
# Save the raster layer
writeRaster(ras, filename = "ras.tif", format="GTiff")
Finally, to use the kriging features from the package gstat
, I need to convert the raster to SpatialPixels
.
# Convert to spatial pixel
st_grid <- rasterToPoints(ras, spatial = TRUE)
gridded(st_grid) <- TRUE
st_grid <- as(st_grid, "SpatialPixels")
st_grid
- this SpatialPixels
that can be used in kriging.
It is an iterative process to determine a suitable grid. Throughout the process, users can change projection, origin, cell size or cell number, depending on their analysis needs.
source to share
If you have an ROI as a polygon imported as SpatialPolygons
, you can use a raster package to rasterize, or use sp::spsample
to sample it using a sample type regular
.
If you do not have such a polygon, you can create points that are regularly spaced along the rectangular length / lat-area, using expand.grid
using seq
lat to generate a sequence of longs and values.
source to share
@yzw and @Edzer make good points for creating a regular rectangular mesh , but sometimes it becomes necessary to create an irregular mesh over a specific polygon, usually for kriging .
This is a rarely documented topic. One good answer can be found here . I am expanding it with the code below:
Consider the built-in dataset meuse. meuse.grid
- irregular mesh. How do we make a grid like meuse.grid for our unique study area?
library(sp) data(meuse.grid) ggplot(data = meuse.grid)+geom_point(aes(x=x, y=y))
Unfortunately, I can't find a built-in spatialPolygonsDataFrame
one to illustrate this next step, so you'll use your own, which I call mySpatialPolygonsDataFrame below. You build a regular rectangular mesh on it, then use the outline of your irregular polygon to anchor this rectangular mesh in an irregular shape.
1. Create a rectangular mesh on top of your SpatialPolygonsDataFrame
grd <- makegrid(mySpatialPolygonsDataFrame, n = 10000)
colnames(grd) <- c('x','y')
2. Get the outline of your irregular polygon
outline <- mySpatialPolygonsDataFrame@polygons[[2]]@Polygons[[1]]@coords
3. Use inout () to clip the rectangular mesh with an outline
require(splancs)
new_grd <- grd[inout(grd,outline), ]
4. Visualize your clipped mesh that can be used for kriging!
ggplot(new_grd) + geom_point(aes(x=x,y=y))
source to share