writeRaster vs write.table或.tif vs .csv

有人給我發送了一些R代碼，使用R包raster等來讀取netCDF文件。該代碼創建一系列*.tif文件。不幸的是，我不是很熟悉raster,*.tif文件或netCDF文件。所以，我試圖改變R的代碼也寫*.csv文件。我認爲下面的代碼以*.tif格式和*.csv格式寫入相同的網格單元格數據。但是，我不確定。我希望有人能夠驗證這兩種格式的數據是相同的。理想情況下，我希望能夠打開*.tif文件並自己進行驗證。我怎樣才能做到這一點？writeRaster vs write.table或.tif vs .csv

如果*.tif文件只包含圖像而不是數字，那麼也許不可能進行直接比較。在這種情況下，我想驗證在*.tif文件中的圖像對應的數據在*.csv文件

下面的R代碼是*.csv文件的內容與有關列和行「後續問題標題。

setwd('c:/users/mark w miller/netCDF/') 

my.file <- "my.netCDF.nc" 
my.var1 <- "my.variable" 

library(ncdf) 
library(rgdal) 
library(chron) 
library(fields) 

file <- open.ncdf(my.file) 
long <- get.var.ncdf(file, varid="lon") 
lat  <- get.var.ncdf(file, varid="lat") 
time <- get.var.ncdf(file, varid="time") 
my.varb <- get.var.ncdf(file, varid=my.var1) 

#netCDF to raster 
library(raster) 

r  <- brick(my.file, varname = my.var1) 

#Crop spatial coverage 
e  <- extent(255,265,35,45) 
rc  <- crop(r, e, bylayer=TRUE) 

lat2 <- lat[ lat >= 35 & lat <= 45] 
long2 <- long[long >= 255 & long <= 265] 

list1 <- unstack(rc) 
rs  <- stack(list1) 

for(i in 1:5){ 
    r2 <- 1+(i-1)*12 
    s2 <- 2+(i-1)*12 
    a2 <- rs[[r2]] 
    b2 <- rs[[s2]] 
    m2 <- stack(a2,b2) 
    my.var <- overlay(m2, fun=function(x,y) {(x+y)}, unstack=TRUE, recycle=FALSE) 
    f2 <- 1999+i 

    writeRaster(my.var, filename=paste("my.var", f2, ".tif"), format="GTiff") 

    my.var2 <- as.matrix(my.var, nrow=length(lat2), byrow=TRUE) 

    write.table(my.var2, file = paste0("my.var", f2, ".csv"), quote = FALSE, sep=",", col.names = FALSE, row.names = FALSE) 
}

這裏有一個*.csv文件的舍入內容：

1.0,0.9,0.8,0.8,0.7,0.7,0.8,0.8,1.0,1.0 
1.0,0.8,0.6,0.5,0.4,0.5,0.7,0.9,1.0,1.0 
1.0,0.7,0.5,0.4,0.3,0.4,0.7,1.0,1.0,1.0 
0.0,0.5,0.4,0.4,0.4,0.6,1.0,1.0,1.0,1.0 
0.0,0.6,0.5,0.4,0.5,0.8,1.0,2.0,2.0,2.0 
1.0,0.7,0.6,0.5,0.6,1.0,1.0,2.0,2.0,2.0 
1.0,0.9,0.8,0.7,0.9,1.0,2.0,2.0,2.0,2.0 
1.0,1.0,1.0,1.0,1.0,2.0,2.0,2.0,1.0,2.0 
2.0,1.0,1.0,1.0,2.0,2.0,2.0,1.0,1.0,2.0 
1.0,1.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,3.0

鑑於：

lat2 
# [1] 44.5 43.5 42.5 41.5 40.5 39.5 38.5 37.5 36.5 35.5 

long2 
# [1] 255.5 256.5 257.5 258.5 259.5 260.5 261.5 262.5 263.5 264.5

我可以放心地下面列和行的名稱分別添加到*.csv文件？

 255.5 256.5 257.5 258.5 259.5 260.5 261.5 262.5 263.5 264.5 
44.5 1.0,0.9,0.8,0.8,0.7,0.7,0.8,0.8,1.0,1.0 
43.5 1.0,0.8,0.6,0.5,0.4,0.5,0.7,0.9,1.0,1.0 
42.5 1.0,0.7,0.5,0.4,0.3,0.4,0.7,1.0,1.0,1.0 
41.5 0.0,0.5,0.4,0.4,0.4,0.6,1.0,1.0,1.0,1.0 
40.5 0.0,0.6,0.5,0.4,0.5,0.8,1.0,2.0,2.0,2.0 
39.5 1.0,0.7,0.6,0.5,0.6,1.0,1.0,2.0,2.0,2.0 
38.5 1.0,0.9,0.8,0.7,0.9,1.0,2.0,2.0,2.0,2.0 
37.5 1.0,1.0,1.0,1.0,1.0,2.0,2.0,2.0,1.0,2.0 
36.5 2.0,1.0,1.0,1.0,2.0,2.0,2.0,1.0,1.0,2.0 
35.5 1.0,1.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,3.0

感謝您的任何建議。實際的NetCDF文件非常大。如果我能弄清楚如何對它進行子集化並將其保存爲相同的netCDF格式，我可能會嘗試將它上傳到某處。

編輯

下面是代碼來創建模擬數據，這些模擬數據轉換爲netCDF文件和分析netCDF文件如在上面的代碼：

setwd('c:/users/mark w miller/netCDF/') 

library(raster) 
library(ncdf) 
library(rgdal) 
library(chron) 
library(fields) 
library(sp) 

set.seed(1234) 

x = seq(255, 269, length = 8) 
y = seq( 36, 40, length = 5) 
xy <- expand.grid(x,y) 

z <- rnorm(nrow(xy), 10, 1) 
rc <- data.frame(xy,z) 

raster.rc1 <- rasterFromXYZ(rc, res=c(2,1), crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) 

z <- rnorm(nrow(xy), 10, 1) 
rc <- data.frame(xy,z) 

raster.rc2 <- rasterFromXYZ(rc, res=c(2,1), crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) 

z <- rnorm(nrow(xy), 10, 1) 
rc <- data.frame(xy,z) 

raster.rc3 <- rasterFromXYZ(rc, res=c(2,1), crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) 

z <- rnorm(nrow(xy), 10, 1) 
rc <- data.frame(xy,z) 

raster.rc4 <- rasterFromXYZ(rc, res=c(2,1), crs=CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")) 


all.years <- list(raster.rc1, raster.rc2, raster.rc3, raster.rc4) 

#all.rc <- stack(all.years) 

all.rc <- brick(all.years) 

writeRaster(all.rc, filename="example.netCDF.nc", format="CDF", bylayer=TRUE, overwrite=TRUE) 


my.file <- open.ncdf('example.netCDF.nc') 
my.file 

long <- get.var.ncdf(my.file, varid="longitude") 
lat <- get.var.ncdf(my.file, varid="latitude") 
time <- get.var.ncdf(my.file, varid="value") 
my.var <- get.var.ncdf(my.file, varid="variable") 

long 
# [1] 255 257 259 261 263 265 267 269 

lat 
# [1] 40 39 38 37 36 

time 
# [1] 1 2 3 4 

my.var 

r  <- brick('example.netCDF.nc', varname = 'variable') 

#Crop spatial coverage 
e  <- extent(257,267,37,39) 
rc  <- crop(r, e, bylayer=TRUE) 

lat2 <- lat[ lat >= 37 & lat <= 39] 
lat2 
long2 <- long[long >= 257 & long <= 267] 
long2 

list1 <- unstack(rc) 
rs  <- stack(list1) 

for(i in 1:2){ 
    r2 <- 1+(i-1)*2 
    s2 <- 2+(i-1)*2 
    a2 <- rs[[r2]] 
    b2 <- rs[[s2]] 
    m2 <- stack(a2,b2) 
    my.sim <- overlay(m2, fun=function(x,y) {(x+y)}, unstack=TRUE, recycle=FALSE) 
    f2 <- 2010+i 

    writeRaster(my.sim, filename=paste("my.sim", f2, ".tif"), format="GTiff") 

    my.sim2 <- as.matrix(my.sim, nrow=length(lat2), byrow=TRUE) 

    write.table(my.sim2, file = paste0("my.sim", f2, ".csv"), quote = FALSE, sep=",", col.names = FALSE, row.names = FALSE) 
}

來源

2014-03-18 Mark Miller

我會從as.data.frame（subset（rc，1），xy = TRUE）開始 – mdsumner

@mdsumner謝謝。這似乎提取了一層數據並將其轉換爲數據幀。也許我可以用'create.ncdf'把這個數據幀轉換成一個netCDF文件。 –

使用writeRaster創建一個netcdf文件 – mdsumner

這是一個普遍的方式比較*.tiff文件和*.csv文件的內容。將示例*.csv文件的內容與其*.tiff情節進行比較，可以讓我確信自己的內容是相同的。

下面我將展示如何顯示*.tiff文件中的數據。

setwd('c:/users/mmiller21/netCDF/') 

library(raster) 

# Here are the contents of 'my.sim2011.csv': 
# 
# 18.31545067 20.22907639 20.34417152 18.11485672 17.93542576 19.52469158 
# 19.20878696 19.43614769 18.41953754 16.42925882 22.05830574 18.31794167 
# 
# compared with the plot of 'my.sim 2011 .tif' 

jpeg(filename = "my.sim.2011.jpeg") 

    r <- raster('my.sim 2011 .tif') 
    plot(r) 
    title(main='my.sim 2011 .tif') 

dev.off()

enter image description here

# Here are the contents of 'my.sim2012.csv': 
# 
# 18.92995739 20.68585968 20.44407845 20.53401566 19.1156435 20.70266819 
# 19.04809856 20.76659107 20.50794601 18.52109146 20.92043018 19.91858768 
# 
# compared with the plot of 'my.sim 2012 .tif' 

jpeg(filename = "my.sim.2012.jpeg") 

    r <- raster('my.sim 2012 .tif') 
    plot(r) 
    title(main='my.sim 2012 .tif') 

dev.off()

enter image description here

這裏是代碼在一個示例*.tiff文件中顯示的數據。該數據與相應的*.csv文件中的數據匹配。

r <- raster('my.sim 2011 .tif') 

r[1,] 
[1] 18.31545 20.22908 20.34417 18.11486 17.93543 19.52469 

r[2,] 
[1] 19.20879 19.43615 18.41954 16.42926 22.05831 18.31794

來源

2014-03-19 05:36:59

writeRaster vs write.table或.tif vs .csv

回答

相關問題