今天接到一位姐妹的求助，她要提取某个站点的PDSI，但是不知道怎么整，问我会不会。这等小事不得让我拿捏一下，我说“坐标，年份，速来”
然后发现用raster包和terra包提取，有一点区别，于是再记录一下。
目的：提取nc格式的站点PDSI
数据源：
Drought indices (uea.ac.uk)
https://crudata.uea.ac.uk/cru/data/drought/
Climate Explorer: Select a monthly field (knmi.nl)
http://climexp.knmi.nl/selectfield_obs2.cgi?id=someone@somewhere
法1：raster包stack处理
library(raster)
nc_stack <- stack("scPDSI.nc")
df <- data.frame(lon=115.983,lat=29.583)
dfsp = SpatialPoints(df)
p <- raster::extract(nc_stack,dfsp)View一下p，发现它长这样
是一个宽数据
class一下p，它是一个矩阵
矩阵、数组
转置一下p
t(p)
转置以后长这样
法2：raster包brick处理
library(raster)
ncc <- brick("scPDSI.nc")
df <- data.frame(lon=115.983,lat=29.583)
dfsp = SpatialPoints(df)
pp <- raster::extract(ncc,dfsp)
t(pp)
brick函数的结果和stack一样
和法1一样，只不过把栅格做成了brick而不是stack，主要区别就在于brick函数秒出结果，stack函数足足等了16多分钟！果然brick函数效率高啊。
法3：terra包rast处理
library(terra)
df <- data.frame(lon=115.983,lat=29.583)
vc <- vect(df)
nc <- rast('scPDSI.nc')
ppp <- terra::extract(nc,vc)
class(ppp)
head(t(ppp))
terra的处理结果
terra的处理速度和raster的brick一样，都是秒出。但是列名（转置后是行名）不一样，rast读取的图层名不如brick函数输出的图层名信息多，brick函数输出的图层名年月日都标上了。这竟然和之前nc转tif的记录中截然相反。
在nc转tif中，rast读取的nc文件图层名是pre_i，brick读取的nc文件名就只是x1，x2等等。
这提了个醒：读取nc文件时要这两种方法都试一试，看看哪个易于后续操作。
法4：还有更快捷的方法，但是只针对于世界气象组织网站上的数据
Climate Explorer: Select a monthly field (knmi.nl)
http://climexp.knmi.nl/selectfield_obs2.cgi?id=someone@somewhere
这个网站收录了很多数据集，包括这次用到的PDSI。
找到需要的数据集
点进去后会要你选择提取的区域，提取点还是面，面的话要加一个mask
经纬度要确定好是十进制的
点击make time series
时间序列图
然后它的时间序列就出来了，还可以下载这个点数据的eps，pdf，元数据，raw data（dat格式），nc文件等等
由于只需要1960到2020年的数据，在这个页面的下面还可以进一步限制条件
好多可以限制的条件
限制完了就可以在新的页面上下载它了。
时间序列
推荐下载dat格式的文件。为什么推荐下载dat格式的文件呢？因为快啊，下载下来的dat格式文件用excel打开，分列，把多余的表头删除就ok了。（别问为啥不用R或者Python或者MATLAB打开dat文件，光是开软件加写路径和代码的工夫excel就在旁边抽事后烟看着你笑了，效率为王啊！）（说到敲文件路径，昨天还和师妹吐槽来着：MATLAB比R人性化太多了，R就是个垃圾，路径还得用反斜杠，右下角的目录窗口更换盘符还得点右边的三个点）
如果不嫌麻烦，就是想用netcdf来install B，其实也可以。用下面的代码：
library(ncdf4)
ls <- nc_open("iscpdsi_115.983E_29.583N_n_1960_2020.nc")
df <- data.frame(ncvar_get(ls))
colnames(df) <- 'PDSI'
library(tidyverse)
year <- rep(1960:2020,each=12)
month <- rep(1:12,times=61)
df <- mutate(df,year=year,month=month)metadata。点这个会下载整个数据集年份的全球区域的原始数据。如果要转区域的tif的话可以参考这个
pdf。下载一个dat文件的pdf版本，属实没必要。
eps。我只在ps上用过eps文件，这方面还真没遇到过，不过图像文件的封装方式有很多，估计这也是其中一种。希望有知道这个文件怎么用遥感方式处理的大手子可以来指点一下。
