使用循环提取坐标,匹配它们并写入文件

问题描述 投票:0回答:2

我试图使用for循环(或apply函数作为替代)从data.frame中提取坐标,搜索E-OBS gridded dataset中的最近点,提取时间x1-x2的温度数据并写入到另一个excel文件。

虽然代码可以提取单个数据点,但我似乎无法在循环中包含此代码,并且可能会在输入坐标旁边添加结果。

library(sp)
library(raster)
library(ncdf4)

#Coordinates
    df
       ID    site                 E        N
1       1   site_place_date1  7.558758 47.81004
2       2   site_place_date2  7.582749 47.63411
3       3   site_place_date3  7.607968 48.01475
4       4   site_place_date4  7.644660 47.67139
       ...     ...   ...              ...`

手动设置目标点的坐标:

lon <- 7.558758  # longitude of location                
lat <- 47.81004 # latitude  of location

#Mean daily temperature
    ncin <- nc_open("tg_0.25deg_reg_v17.0.nc")
      print(ncin)
      t <- ncvar_get(ncin,"time")
      tunits <- ncatt_get(ncin,"time","units")nt <- dim(t)
      nt
      obsoutput <- ncvar_get(ncin, 
                       start= c(which.min(abs(ncin$dim$longitude$vals -   lon)), # look for closest long
                                which.min(abs(ncin$dim$latitude$vals -  lat)),  # look for closest lat
                                1),
                       count=c(1,1,-1))
      DataMeanT <- data.frame(DateN= t, MeanDailyT = obsoutput)
      nc_close(ncin)
      head(DataMeanT)


#check if there are NAs =999
    summary(DataMeanT)

    Data = DataMeanT
    Data$Date = as.Date(Data$DateN,origin="20000-01-01")
    Data$Year = format(Data$Date,"%Y")
    Data$Month = format(Data$Date,"%m")
    head(Data)
    Data$YearMonth = format(Data$Date, format="%Y-%b")

    Data_annual = aggregate(("T_AnnualMean" = MeanDailyT) ~ Year,data = Data, FUN = mean,na.action = na.pass)
    names(Data_annual)[2] <- "AirT"
    head(Data_annual)

#Export table
    write.table(Data_annual, "Site_AirTemp.csv", row.names = FALSE, append = FALSE, col.names = TRUE, sep = ", ", quote = TRUE)

目的是将脚本作为df中所有坐标的循环的一部分运行,并将温度数据写入具有site-ID信息的新数据表,或者写入df的下一列。

r loops for-loop apply sapply
2个回答
0
投票

只需将整个过程包装在一个已定义的方法中,然后使用apply函数传入lon / lat坐标。一个伟大的候选人是mapply或其包装Mapdf$Edf$N的两个向量之间迭代元素。此外,第三个参数df$site被传递给唯一CSV名称的方法,因为现在同一个文件将被覆盖。

下面的一些非赋值行如headsummary被删除,因为它们在方法中什么都不做。此外,上下文管理器withinwith用于避免重复使用Data$以实现更精简的数据操作。 Map调用写入文件并构建聚合数据帧列表供以后使用。

功能

my_function <- function(lon, lat, site) {    
    # Mean daily temperature
    ncin <- nc_open("tg_0.25deg_reg_v17.0.nc")
      print(ncin)
      t <- ncvar_get(ncin,"time")
      tunits <- ncatt_get(ncin,"time","units")nt <- dim(t)

      # look for closest lon and lat
      obsoutput <- ncvar_get(ncin, 
                             start = c(which.min(abs(ncin$dim$longitude$vals - lon)),
                                      which.min(abs(ncin$dim$latitude$vals - lat)),
                                      1),
                             count = c(1,1,-1))

      DataMeanT <- data.frame(DateN = t, MeanDailyT = obsoutput)
    nc_close(ncin)    

    Data <- within(DataMeanT, {
               Date <- as.Date(DateN, origin="2000-01-01")
               Year <- format(Date,"%Y")
               Month <- format(Date,"%m")
               YearMonth <- format(Date, format="%Y-%b")
            })

    Data_annual <- with(Data, aggregate(list("AirT" = MeanDailyT), list(Year=Year),
                                        FUN = mean, na.action = na.pass))    
    # Export table
    write.table(Data_annual, paste0("Site_AirTemp_", site, "_.csv"), row.names=FALSE,
                append = FALSE, col.names = TRUE, sep = ", ", quote = TRUE)

    # SAVE AGGREGATED DATA FRAME
    return(Data_annual)
}

呼叫

# ITERATE THROUGH EACH LON/LAT PAIR ELEMENTWISE
df_list <- Map(my_function, df$E, df$N, df$site)

# df_list <- mapply(my_function, df$E, df$N, df$site, SIMPLIFY=FALSE)    # EQUIVALENT CALL

0
投票

很难回答你的问题,因为它不可复制。但你可能会这样做:

library(raster)
b <- brick("tg_0.25deg_reg_v17.0.nc")
e <- extract(b, df[, c('E', 'N')])
© www.soinside.com 2019 - 2024. All rights reserved.