地理空间数据一般都带有坐标系,最常用的是WGS 84,如果我们已经读入了一个地理数据,比如保存在data中,那么可以查看它的空间坐标系:
library(pacman)
p_load(sf)
data = st_read("G:/Rdata/China/fujian1.shp")
# Reading layer `fujian1' from data source `G:\Rdata\China\fujian1.shp' using driver `ESRI Shapefile'
# Simple feature collection with 10 features and 2 fields
# geometry type: POLYGON
# dimension: XY
# bbox: xmin: 1115192 ymin: 3088867 xmax: 1541805 ymax: 3636025
# projected CRS: China_Lambert_Conformal_Conic
st_crs(data)$proj4string
# [1] "+proj=lcc +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs"
如果数据中没有坐标系,那么我们可以通过st_set_crs函数来指定它的坐标系,坐标系的表示方法有两种:
1、Proj4:由一长串字符串构成,比如"+proj=lcc +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs"
2、EPSG:由数字编码构成,如4326(WGS 84)后者要更方便一些,比如对于一个没有坐标系的数据,我们可以这样设置其坐标系以WGS84坐标系(4326)为例。
data <- data %>%
st_set_crs(., 4326)
# Warning message:
# st_crs<- : replacing crs does not reproject data; use st_transform for that
#因为有坐标,所以提示使用st_transform以转到WGS84坐标系。
如果数据有坐标系,需要转换坐标系的话,可以用st_transfrom函数
dataProjected <- data %>%
st_transform(.,4326)
st_crs(dataProjected)$proj4string
> st_crs(data)$proj4string
[1] "+proj=lcc +lat_0=0 +lon_0=105 +lat_1=30 +lat_2=62 +x_0=0 +y_0=0 +ellps=krass +units=m +no_defs"
> st_crs(dataProjected)$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"
通常我们读取数据为sf格式,sf和sp数据可以通过一定的格式转换,例子如下,假设Kenoutline为我们的地理数据:
#From sf to sp
KenoutlineSP <- Kenoutline %>%
as(., "Spatial")
#From sp to sf
KenoutlineSF <- KenoutlineSP %>%
st_as_sf()
对于栅格数据,则需要使用raster包来加载,raster()函数可以载入数据,stack()可以做波段叠加,而projectRaster()函数则可以根据其proj4的字符串来进行不同形式的投影,
p_load(raster)
Aug <- raster("G:/Rdata/wc2.1_10m_tavg/wc2.1_10m_tavg_08.tif")
Aug
# class : RasterLayer
# dimensions : 1080, 2160, 2332800 (nrow, ncol, ncell)
# resolution : 0.1666667, 0.1666667 (x, y)
# extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
# crs : +proj=longlat +datum=WGS84 +no_defs
# source : G:/Rdata/wc2.1_10m_tavg/wc2.1_10m_tavg_08.tif
# names : wc2.1_10m_tavg_08
# values : -66.5225, 38.43275 (min, max)
plot(Aug)
plot(Aug)
newProj <- "+proj=leac"
# get the Aug raster and give it the new proj4
pro1 <- Aug %>%
projectRaster(.,crs=newProj)
plot(pro1)
newProj <- "+proj=lcca +lat_0=35"
# get the Aug raster and give it the new proj4
pro1 <- Aug %>%
projectRaster(.,crs=newProj)
plot(pro1)
newProj <- "+proj=aeqd"
# get the Aug raster and give it the new proj4
pro1 <- Aug %>%
projectRaster(.,crs=newProj)
plot(pro1)
newProj <- "+proj=aea +lat_1=29.5 +lat_2=42.5"
# get the Aug raster and give it the new proj4
pro1 <- Aug %>%
projectRaster(.,crs=newProj)
plot(pro1)
plot(Aug)
newProj <- "+proj=gn_sinu +m=2 +n=3"
# get the Aug raster and give it the new proj4
pro1 <- Aug %>%
projectRaster(.,crs=newProj)
plot(pro1)