R语言空间数据分析学习笔记3——坐标系设置

地理空间数据一般都带有坐标系，最常用的是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)