How to georeference images using R
Juan Carlos Rubio-Polania
2023-12-02
Georeferencing involves assigning geographic coordinates, typically latitude and longitude, to data, anchoring it to specific locations on the Earth’s surface. This practice is essential in diverse fields, including geography, environmental science, urban planning, and agriculture, enabling accurate representation on maps and enhancing spatial analysis.
Georeferencing images is a crucial aspect, serving as a foundation for mapping and integrating spatial data. It allows for the alignment of image features with real-world locations, providing spatial context that enhances data understanding. Notably, in applications like remote sensing, georeferenced images unlock valuable information not previously available, playing a pivotal role in monitoring environmental changes and assessing land cover. The precision of georeferencing supports comprehensive spatial analysis and contributes to navigation and positioning systems, facilitating GPS-based location tracking and precise object identification on the Earth’s surface.
Presented here is the method for georeferencing images using R and Python, based on a map provided by Mendoza et al. in 2023. The starting point of this project is a screenshot of Figure 1 from the scientific paper, which discusses research conducted at the lakes across the Pyrenees.
Later, to capture the screenshot, it is loaded into R using the raster function of the raster package
IMA = raster('C:/Pyrenees.png')Visualize the map using rastervis package
levelplot(IMA,col.regions = viridis::viridis(72))
Fig. 1: Print screen figure loaded using the raster package and plotted using rasterVis packages
Correcting the image using strech function of raster package
The stretch function in the R raster package is used to adjust or change the scale of values in a raster image. This function redistributes pixel values in a raster image to enhance visualization by highlighting specific value ranges and improving contrast. You can use the parameters min_value and max_value, or use the predefined ones in the package. The figures are visualized using the levelplot function of the rasterVis package
image2 = stretch(IMA)
levelplot(image2,col.regions = viridis::viridis(72))Fig 2. Print screen figure adjusted using the stretch function
Defining the extension using extent function of raster package
The extent function in the R raster package is used to obtain or modify the spatial extent of a raster object, defining the geographical region covered by the raster. You can use the code below to retrieve the current extent.
To set a new extent using the minimum and maximum coordinates on the x and y axes, you can utilize the code presented below.
In the case of our example:
extent(image2) = c(-2.01, 3.25, 41.8, 43.25)Adding projection
The projection function in the R raster package is used to obtain or modify the spatial projection of a raster object, defining how spatial data, such as coordinates, relates to the Earth’s surface. This function is useful when working with raster data in specific geospatial contexts, as different analyses and visualizations may require particular projections. You can use to obtain the current projection:
To set a new projection, use ‘NewProj’ as the code or definition for the projection you want to assign to the raster object.
The image’s projection is specified here using a common projection for mountains and high latitudes (Lambert).
projection(image2) = crs("+proj=lcc +lat_1=41.8 +lat_2=43.25 +lat_0=42.525 +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs") To visualize the image, layers were stacked, and then a data frame was created. Subsequently, its column names were changed to construct a figure using the RGB format and ggplot functions.
image3 = stack(image2,image2,image2)
DF = as.data.frame(image3, xy = TRUE)
head(DF)## x y layer.1 layer.2 layer.3
## 1 -2.007531 43.2484 84 84 84
## 2 -2.002592 43.2484 19 19 19
## 3 -1.997653 43.2484 21 21 21
## 4 -1.992714 43.2484 19 19 19
## 5 -1.987775 43.2484 19 19 19
## 6 -1.982836 43.2484 19 19 19colnames(DF) = c('x', 'y', 'Red', 'Green', 'Blue')
head(DF)## x y Red Green Blue
## 1 -2.007531 43.2484 84 84 84
## 2 -2.002592 43.2484 19 19 19
## 3 -1.997653 43.2484 21 21 21
## 4 -1.992714 43.2484 19 19 19
## 5 -1.987775 43.2484 19 19 19
## 6 -1.982836 43.2484 19 19 19ggplot(data = DF, aes(x = x, y =y))+ #plot map
geom_raster(fill = rgb(DF$Red/255, DF$Green/255, DF$Blue/255)) +
scale_y_continuous(expand=c(0,0)) + scale_x_continuous(expand=c(0,0)) +
theme(text= element_text(size=25, family = "Arial", colour = "black")) +
labs(x = 'Longitude', y = 'Latitude') +
theme_classic()
Fig 3. Georeferended data got from locator function
Despite the appearance of the figure, the map maintains reliable georeferencing
raster_df = as.data.frame(image2, xy = TRUE)
world = map_data("world")Finally, create the georeferenced file using writeRaster function
writeRaster(image2,'C:/Users/jucar/Documents/GitHub/Rpubs/PyreneesRF.tif',
format = "GTiff", overwrite = TRUE)Load the georeferenced file into R
The georeferenced file is loaded into R, and the stacked file is created. Additionally, a new data frame is constructed following the RGB format mentioned above. Later, the figures were plotted using the ggplot2 package.
img = brick('C:/PyreneesRF.tif')
img## class : RasterBrick
## dimensions : 453, 1065, 482445, 1 (nrow, ncol, ncell, nlayers)
## resolution : 0.004938967, 0.003200883 (x, y)
## extent : -2.01, 3.25, 41.8, 43.25 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=42.525 +lon_0=0 +lat_1=41.8 +lat_2=43.25 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
## source : PyreneesRF.tif
## names : PyreneesRF
## min values : 0
## max values : 255img2 = stack(img,img,img)
img2## class : RasterStack
## dimensions : 453, 1065, 482445, 3 (nrow, ncol, ncell, nlayers)
## resolution : 0.004938967, 0.003200883 (x, y)
## extent : -2.01, 3.25, 41.8, 43.25 (xmin, xmax, ymin, ymax)
## crs : +proj=lcc +lat_0=42.525 +lon_0=0 +lat_1=41.8 +lat_2=43.25 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs
## names : PyreneesRF.1, PyreneesRF.2, PyreneesRF.3
## min values : 0, 0, 0
## max values : 255, 255, 255DF2 = as.data.frame(img2, xy = TRUE)
head(DF2)## x y PyreneesRF.1 PyreneesRF.2 PyreneesRF.3
## 1 -2.007531 43.2484 84 84 84
## 2 -2.002592 43.2484 19 19 19
## 3 -1.997653 43.2484 21 21 21
## 4 -1.992714 43.2484 19 19 19
## 5 -1.987775 43.2484 19 19 19
## 6 -1.982836 43.2484 19 19 19colnames(DF2) = c('x', 'y', 'Red', 'Green', 'Blue')
head(DF2)## x y Red Green Blue
## 1 -2.007531 43.2484 84 84 84
## 2 -2.002592 43.2484 19 19 19
## 3 -1.997653 43.2484 21 21 21
## 4 -1.992714 43.2484 19 19 19
## 5 -1.987775 43.2484 19 19 19
## 6 -1.982836 43.2484 19 19 19ggplot(data = DF2, aes(x = x, y =y))+ #plot map
geom_raster(fill = rgb(DF2$Red/255, DF2$Green/255, DF2$Blue/255)) +
scale_y_continuous(expand=c(0,0)) + scale_x_continuous(expand=c(0,0)) +
theme(text= element_text(size=25, family = "Arial", colour = "black")) +
labs(x = 'Longitude', y = 'Latitude') +
theme_bw()
Fig. 4: Georeferenced data obtained from the locator function
Capturing the location of each station using locator function
The positions are obtained by clicking with the mouse on the location represented or marked on the map.
Getting the information obtained from locator
In the R Markdown document, the locator function is causing issues. The data presented here was obtained using the function in an R script.
df1 = data.frame(x1,y1)
colnames(df1) = c('Longitude', 'Latitude')
head(df1)## Longitude Latitude
## 1 -0.6311033 42.75157
## 2 -0.4157956 42.81119
## 3 -0.4124832 42.71182
## 4 -0.2832986 42.74826
## 5 -0.2071128 42.57932
## 6 -0.2104252 42.64888Now, it is possible to plot the selected locations.
Fig. 5: Georeferenced data obtained from the locator function
The result obtained will be to depend of the how reliable are the extreme values used in the extent function and the projection used that in this case was Lambert projection because this area is located between 41.8°N and 43.5°N and is a highland area.
References
De Mendoza, G., Araujo, R., & Catalan, J. (2023). Factors Influencing the Distribution of Freshwater Mollusks in the Lakes of the Pyrenees: Implications in a Shifting Climate Scenario. Diversity, 15(4), 500. https://doi.org/10.3390/d15040500
Deckmyn, A. (2018). Draw Geographical Maps. Package ‘maps’ [R language].
Delignette-Muller, M. L., Pouillot, R., & Denis, J.-B. (2008). Use of the library fitdistrplus to specify a distribution from data [R language].
Hijmans, R., & Elith, J. (2011). Species distribution modelling with R (p. 87).
Hijmans, R. J. (2016). Writing functions with the ”raster” package (p. 11).
Hijmans, R. J. (2019). Introduction to the ”geosphere” package [R language].
Kahle, D., & Wickham, H. (2013). ggmap: Spatial Visualization with ggplot2. The R Journal, 5(1), 144. https://doi.org/10.32614/RJ-2013-014
Ushey, K., Allaire, J., Tang, Y., Eddelbuettel, E., Lewis, B., Keydana, S., Hafen, R., & Geelnard, M. (2020). Interface to “Python”. Package ‘reticulate.’
Zaniewski, A. E., Lehmann, A., & Overton, J. M. (2002). Predicting species spatial distributions using presence-only data: A case study of native New Zealand ferns. Ecological Modelling, 157(2–3), 261–280. https://doi.org/10.1016/S0304-3800(02)00199-0