将 ASCII 格式的海面温度文本文件导入 R

Posted 2023-04-18

【中文标题】将 ASCII 格式的海面温度文本文件导入 R

【英文标题】：Importing Sea Surface Temperature text files in ASCII format into R

【发布时间】：2014-06-09 06:54:13

【问题描述】：

我已经为Hadley Sea Surface Temperature observations 下载了多个 .txt.gz 文件。数据已解压缩，导致ASCII format 中有多个.txt 文件。 我有以下文件（R 脚本是我正在处理的文件）：

list.files()
 [1] "Get_SST_Data.R"                "HadISST1_SST_1931-1960.txt"    "HadISST1_SST_1931-1960.txt.gz"
 [4] "HadISST1_SST_1961-1990.txt"    "HadISST1_SST_1961-1990.txt.gz" "HadISST1_SST_1991-2003.txt"   
 [7] "HadISST1_SST_2004.txt"         "HadISST1_SST_2005.txt"         "HadISST1_SST_2006.txt"        
[10] "HadISST1_SST_2007.txt"         "HadISST1_SST_2008.txt"         "HadISST1_SST_2009.txt"        
[13] "HadISST1_SST_2010.txt"         "HadISST1_SST_2011.txt"         "HadISST1_SST_2012.txt"        
[16] "HadISST1_SST_2013.txt"    

我希望能够利用温度数据为自 1950 年以来每天的海面温度制作一个数字向量，最终制作一个时间序列图。

看起来像这样

[附注]仅供参考...]

提前致谢！

【问题讨论】：

看起来像this。

我在您的问题本身中添加了您在上述评论中提到的图像...

谢谢@vrajs5，我做不出来

顺便说一句，您的问题是创建图表还是导入多个文件？无论如何，类似的问题已经得到解答......

不，最终结果是在模型中使用温度变量。我只是不确定如何将数据组织成有意义的格式。 ASCII 数据似乎不可读。我已经导入了文件，没问题

【参考方案1】：

NetCDF 绝对是一个更好的方法，因为 ascii 数据的格式非常糟糕。也就是说，这是一个读取您下载的数据的函数。

read.things <- function(f) 
  # f is the file path of your ascii data
  require(raster)
  d <- readLines(f)
  d <- split(d, rep(1:12, each=181))
  d <- lapply(d, function(x) read.fwf(textConnection(x), rep(6, 360), 
                                      skip=1, stringsAsFactors=FALSE,
                                      na.strings=c(-1000, -32768)))
  d <- lapply(d, function(x) sapply(x, as.numeric))
  out <- stack(lapply(d, raster))
  names(out) <- month.abb
  extent(out) <- c(-180, 180, -90, 90)
  out/100

请注意，我已将 100% 冰格 (-100) 和陆地格 (-32768) 设置为 NA。

下面，我们以下载其中一个文件（1Mb）为例：

download.file(
  'http://www.metoffice.gov.uk/hadobs/hadisst/data/HadISST1_SST_2004.txt.gz',
  destfile= f <- tempfile())

s <- read.things(f)

s

# class       : RasterBrick 
# dimensions  : 180, 360, 64800, 12  (nrow, ncol, ncell, nlayers)
# resolution  : 1, 1  (x, y)
# extent      : -180, 180, -90, 90  (xmin, xmax, ymin, ymax)
# coord. ref. : NA 
# data source : in memory
# names       :   Jan,   Feb,   Mar,   Apr,   May,   Jun,   Jul,   Aug,   Sep,   Oct,   Nov,   Dec 
# min values  :   -10,   -10,   -10,   -10,   -10,   -10,   -10,   -10,   -10,   -10,   -10,   -10 
# max values  : 30.34, 30.58, 30.43, 30.50, 30.83, 31.39, 32.71, 33.40, 32.61, 31.52, 30.60, 30.51 

library(rasterVis)
levelplot(s, at=seq(min(s[], na.rm=T), max(s[], na.rm=T), len=100),
          col.regions=colorRampPalette(c('#2c7bb6', '#abd9e9', '#ffffbf', 
                                         '#fdae61', '#d7191c')))

【讨论】：

这太可爱了！但是，我应该澄清一下 - 我只想要过去 60 年特定区域（纬度 = 35，经度 = 116）的温度数据。期望的结果类似于与温度对齐的时间/日期变量。我正在使用海面温度与时间序列进行比较。虽然情节很美

@eboylen 我假设您的意思是 lat = -35，否则您将在中国的陆地牢房中（我注意到您来自西澳）。不管怎样，一旦你得到了rasterStack，你就可以这样做：s[cellFromXY(s, c(116, -35))]，它返回包含该点的网格单元的每月 SST。

感谢@jbaums，这适用于只有 1 年数据的数据集。但是，有三个数据集包含 30 年的数据，例如 HadISST1_SST_1931-1960.txt.gz - 有什么建议如何处理这些数据集？

【参考方案2】：

[编辑：仅在 Linux 上测试]

R 能够读取 NetCDF 格式 (http://www.metoffice.gov.uk/hadobs/hadisst/data/HadISST_sst.nc.gz)。解压后可以使用“raster”包读取这些数据，如：

library(raster)
library(xts)
library(caTools)    

一些时间定义：

startYear <- 1950   # start of the period
endYear <- 2011     # end of the period
subp <- '1951-01-01/1980-12-01'   # period for the climatology calculation

打开文件：

sst <- brick('HadISST_sst.nc')
Date <- substr(names(sst),2,11) 
Date <- gsub('\\.', '\\-', Date)
Date <- as.Date(Date)
dstart <- paste(startYear,'01','01',sep='-'); dstart <- grep(dstart, Date)
dend <- paste(endYear,'12','01',sep='-'); dend <- grep(dend, Date)
sst <- subset(sst, dstart:dend)
Date <- Date[dstart:dend]

提取特定点的时间序列（lat=35，lon=120）：

tserie <- as.vector(extract(sst, cbind(116, -35)))
tserie <- xts(tserie, order.by=Date)

计算 subp 期的气候学：

clim <- as.numeric()
for(ii in 1:12)
  clim[ii] <- mean(tserie[subp][(.indexmon(tserie[subp])+1) == ii])

clim <- xts(rep(clim, length(tserie)/12), order.by=Date)

计算异常：

tserie <- tserie - clim

绘制结果：

par(las=1)
plot(tserie, t='n', main='HadISST')
lines(tserie, col='grey')
lines(xts(runmean(tserie, 12), order.by=Date), col='red', lwd=2)
legend('bottomleft', c('Monthly anomaly','12-month moving avg'), lty=c(1,1), lwd=c(1,2), col=c('grey','red'))

【讨论】：

我对 sst 矩阵变量中的内容有点困惑

这看起来很棒@Pascal - 但是使用dstart <- grep(dstart, Date) dend <- grep(dend, Date) 命令我遇到了一些麻烦。他们似乎在制作空整数 - 这导致运行 sst <- subset(sst, dstart:dend) 时出现错误 sst <- subset(sst, dstart:dend)

我没有收到此错误。你设置startYear和endYear了吗？

是的，我做到了。发生的情况是第一个 dstart / dend 设置日期，然后 grep 将其覆盖为新的空数字变量

它对我来说很好用。在这个例子中，dstart 给出了961，dend 给出了1704。

【参考方案3】：

您会收到错误消息，因为如果您查看日期的结构，它们总是从 16 日到 16 日。如果你更换：

dstart <- paste(startYear,'01','01',sep='-'); dstart <- grep(dstart, Date)
dend <- paste(endYear,'12','01',sep='-'); dend <- grep(dend, Date)

为，

dstart <- paste(startYear,'01','16',sep='-'); dstart <- grep(dstart, Date)
dend <- paste(endYear,'12','16',sep='-'); dend <- grep(dend, Date)

它会起作用的。