根据 .txt 中的纬度/经度提取时间序列

Question

之前已经回答了与此类似的问题，但我找不到解决我的特定问题的方法：

我有一个包含辐照度值的 .txt，我想在 csv 文件中提取特定纬度和经度的辐照度值的时间序列。 .txt 文件的（头）看起来像：

NCOLS 839
NROWS 679
XLLCORNER 112.025
YLLCORNER -43.975
CELLSIZE 0.05
NODATA_VALUE -999
-999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999 -999  -999 -999 -999 -999...

该文件没有 lon 和 lat 列（这是提取数据的直接方法）。然而它有一个像这样的标题：NCOLS，NROWS，XLLCORNER，......但我不知道如何提取一个看起来像这样的文件的经纬度系列。

此外，一旦提取了特定纬度/经度的时间序列，我想将可能的 nodata (-999) 替换为该系列的前一个值和下一个值的平均值。

我使用的是 R 版本 3.2.4。

指向 .txt 文件之一的可下载链接。请注意此文件适用于某个年/月/日/小时。一年中有几个小时的文件。 http://www.mediafire.com/download/t7xl9adxz3764f3/solar_dni_20130102_00UT.txt

Answer 1

该文件是一个ESRI ASCII Raster format 文件。将其名称更改为 solar.asc（.asc 似乎 - 不幸的是 - 很重要），然后安装 SDMTools 包。完成后，您可以：

library(SDMTools)

dat <- read.asc("solar.asc")

或

df <- asc2dataframe("solar.asc")

绝对是栅格数据（不是时间序列）。

更新

因为您需要一个完整的演练。

library(SDMTools)

# make sure you're in the same directory as your file
getwd()

## [1] "/Path/to/the/directory/with/the/data"

# make sure your original file is in the directory
file.exists("solar_dni_20130102_00UT.txt")

## [1] TRUE

# since the SDMTools functions seem to need ".asc" at the end
file.rename("solar_dni_20130102_00UT.txt",
            "solar_dni_20130102_00UT.asc")

# this reads it as a raster
dat <- read.asc("solar_dni_20130102_00UT.asc")

# check it out
str(dat)

##  asc [1:839, 1:679] NA NA NA NA NA NA NA NA NA NA ...
##  - attr(*, "xll")= num 112
##  - attr(*, "yll")= num -44
##  - attr(*, "cellsize")= num 0.05
##  - attr(*, "type")= chr "numeric"


# this reads it as a data.frame
df <- asc2dataframe("solar_dni_20130102_00UT.asc")

## working with solar_dni_20130102_00UT.asc ...

# check it out
str(dat)

##  asc [1:839, 1:679] NA NA NA NA NA NA NA NA NA NA ...
##  - attr(*, "xll")= num 112
##  - attr(*, "yll")= num -44
##  - attr(*, "cellsize")= num 0.05
##  - attr(*, "type")= chr "numeric"

# put things back the way they were
file.rename("solar_dni_20130102_00UT.asc",
            "solar_dni_20130102_00UT.txt")

而且，我们甚至可以证明这是一个 ESRI ASCII 光栅：

quick.map(dat, 0.5)

哇！澳大利亚（文件中的纬度/经度已经泄露）。

更新 #2

您可以通过以下方式简化提取：

library(raster)

r <- raster.from.asc(dat)

as.numeric(extract(r, cbind(150.757, -26.918)))
## [1] 973