今天接到一位姐妹的求助，她要提取某個站點的PDSI，但是不知道怎么整，問我會不會。這等小事不得讓我拿捏一下，我說“坐標，年份，速來”

然后發(fā)現(xiàn)用raster包和terra包提取，有一點區(qū)別，于是再記錄一下。

目的：提取nc格式的站點PDSI

數(shù)據(jù)源：

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處理

View一下p，發(fā)現(xiàn)它長這樣

class一下p，它是一個矩陣

轉(zhuǎn)置一下p

法2：raster包brick處理

和法1一樣，只不過把柵格做成了brick而不是stack，主要區(qū)別就在于brick函數(shù)秒出結(jié)果，stack函數(shù)足足等了16多分鐘！果然brick函數(shù)效率高啊。

法3：terra包rast處理

terra的處理速度和raster的brick一樣，都是秒出。但是列名（轉(zhuǎn)置后是行名）不一樣，rast讀取的圖層名不如brick函數(shù)輸出的圖層名信息多，brick函數(shù)輸出的圖層名年月日都標上了。這竟然和之前nc轉(zhuǎn)tif的記錄中截然相反。

在nc轉(zhuǎn)tif中，rast讀取的nc文件圖層名是pre_i，brick讀取的nc文件名就只是x1，x2等等。

這提了個醒：讀取nc文件時要這兩種方法都試一試，看看哪個易于后續(xù)操作。

法4：還有更快捷的方法，但是只針對于世界氣象組織網(wǎng)站上的數(shù)據(jù)

Climate Explorer: Select a monthly field (knmi.nl)

http://climexp.knmi.nl/selectfield_obs2.cgi?id=someone@somewhere

這個網(wǎng)站收錄了很多數(shù)據(jù)集，包括這次用到的PDSI。

點進去后會要你選擇提取的區(qū)域，提取點還是面，面的話要加一個mask

點擊make time series

然后它的時間序列就出來了，還可以下載這個點數(shù)據(jù)的eps，pdf，元數(shù)據(jù)，raw data（dat格式），nc文件等等

由于只需要1960到2020年的數(shù)據(jù)，在這個頁面的下面還可以進一步限制條件

限制完了就可以在新的頁面上下載它了。

• 推薦下載dat格式的文件。為什么推薦下載dat格式的文件呢？因為快啊，下載下來的dat格式文件用excel打開，分列，把多余的表頭刪除就ok了。（別問為啥不用R或者Python或者MATLAB打開dat文件，光是開軟件加寫路徑和代碼的工夫excel就在旁邊抽事后煙看著你笑了，效率為王?。。ㄕf到敲文件路徑，昨天還和師妹吐槽來著：MATLAB比R人性化太多了，R就是個垃圾，路徑還得用反斜杠，右下角的目錄窗口更換盤符還得點右邊的三個點）

• 如果不嫌麻煩，就是想用netcdf來install B，其實也可以。用下面的代碼：

• metadata。點這個會下載整個數(shù)據(jù)集年份的全球區(qū)域的原始數(shù)據(jù)。如果要轉(zhuǎn)區(qū)域的tif的話可以參考這個

• pdf。下載一個dat文件的pdf版本，屬實沒必要。

• eps。我只在ps上用過eps文件，這方面還真沒遇到過，不過圖像文件的封裝方式有很多，估計這也是其中一種。希望有知道這個文件怎么用遙感方式處理的大手子可以來指點一下。