1

在过去的三周里，我只使用 R，所以我还是个新手。我目前正在使用来自新英格兰（路径）的 NetCDF 气候数据。我还有一个 .csv 文件，其中包含我们要查看的特定城市的坐标（cities.path）。我已经能够从 .csv 文件中与特定城市对应的模型网格单元中提取时间年度时间序列和趋势。问题在于能够将我的 .csv 文件中的重要城市车站以年平均值绘制到地图上。

当我运行脚本行 'ave_annual_cities <- extract(annual_ave_stack, city.points, df = T)' 时，我得到一个包含重要城市点的纬度/经度图。当我在控制台中运行 'plot(coordinates(cities.points))' 时，我再次得到一个包含重要城市点的纬度/经度图，但它是独立的，与我的 ave_annual_cities 图分开。当我运行 'levelplot(subset(annual_ave_stack, 1), margin=F) + layer(sp.points(cities.points, plot.grid = TRUE))' 时，我得到一张新英格兰的年平均值图表。

到目前为止，这是我的脚本。

#Coordinate CS lat/lon pull from Important City Stations
# imports the csv of lat/lon for specific cities
# reads the lat/lon of the .nc modeled climate files
# extracts the time annual time series and trend from model grid cells 
corresponding to your specific cities' locations.
#Graphs City points according to annual time series graph

# Libraries (probably do not need all)
library(survival)
library(lattice)
library(ncdf4)
library(RNetCDF)
library(date)
library(raster)
library(rasterVis)
library(latticeExtra)
library(RColorBrewer)
library(rgdal)
library(rgeos)
library(wux)

path <- "/net/nfs/merrimack/raid/Northeast_US_Downscaling_cmip5/"
vars = c("tasmin","tasmax")
mods = c("ACCESS1-0", "ACCESS1-3",
     "bcc-csm1-1-m", "bcc-csm1-1")
scns = c("rcp45", "rcp85")

cities.path <-
"/net/home/cv/marina/Summer2017_Projects/Lat_Lon/NE_Important_Cities.csv"
necity.vars <- c("City", "State", "Population", "Latitude", "Longitude", 
"Elevation(meters")

#These both read the .csv file (first uses 'utils', second uses 'wux')
#1
cities.read <- read.delim(cities.path, header = T, sep = ",")
#2
read.table <- read.wux.table(cities.path)
cities.read <- subset(cities.read, subreg = "City", sep = ",")

# To test one coordinate point
point_1 <- c("test.city", 44.31, -69.79)
colnames(point_1)<-c("cities", "latitude", "longitude" )

# To read only "Cities", "Latitude", and "Longitude"
cities.points <- subset(cities.read, select = c(1, 4, 5))
cities.points <- as.data.frame(cities.points)
colnames(cities.points)<- c("City", "Latitude", "Longitude" )

#Set plot coordinates for .csv graph
coordinates(cities.points) <- ~ Longitude + Latitude
proj4string(cities.points) <- c("+proj=longlat +datum=WGS84 +ellps=WGS84 
+towgs84=0,0,0") 

# Start loop to envelope all .nc files
for (iv in 1:2){
  for (im in 1:4){
    for (is in 1:2){
      for(i in 2006:2007){
        full<-paste(path, vars[iv], "_day_", mods[im], "_", scns[is], 
"_r1i1p1_", i, "0101-", i, "1231.16th.nc", sep="")
    # this will print out 
    #/net/nfs/merrimack/raid/Northeast_US_Downscaling_cmip5/NameOfFiles.nc

    # this line will print the full file name
    print(full)

    # use the brick function to read the full netCDF file.
    # note: the varname argument is not necessary, but if a file has multiple 
varables, brick will read the first one by default.
    air_t<-brick(full, varname=vars[iv])

    # use the calc function to get an average for each grid cell over the 
entire year
    annual_ave_t<-calc(air_t, fun = mean, na.rm=T)

        if(i == 2006){
          annual_ave_stack = annual_ave_t
        }else{
          annual_ave_stack<-stack(annual_ave_stack, annual_ave_t)
        }  # end of if/else 
      }   # end of year loop
      #extract annual means for grid cells for each year corresponding to 
 important cities
  ave_annual_cities <- extract(annual_ave_stack, cities.points, df = T)
    }   # end of scenario loop
  }   # end of model loop
}  # end of variable loop


levelplot(subset(annual_ave_stack, 1), margin=F) + 
  layer(sp.points(cities.points, plot.grid = TRUE))

# Read lat/lon from .nc climate files
# http://geog.uoregon.edu/bartlein/courses/geog607/Rmd/netCDF_01.htm
climate.data <- nc_open(full)
lat <- ncvar_get(climate.data, varid = "lat")
nlat <- dim("lat")
lat
lon <- ncvar_get(climate.data, varid =  "lon")
nlon <- dim("lon")
lon
  # This gives all lat data. 

#print long and lat variables: confirms the dimensions of the data
print(c(nlon, nlat))


  # If I need time series... 
my.time <- nc_open(climate.data, "time")
n.dates <- trunc(length(climate.data))
n.dates

# open NetCDF choosing pop and lat/lon points
cities.pop.points <- subset(cities.read, select = c(1, 3, 4, 5))
# print NetCDF coordinates with pop
print(cities.pop.points)

希望这是有道理的。

4

1 回答 1

0

如果您的输入数据被正确读取为多波段栅格，您可以使用以下命令将 netcdf 作为栅格读取：

 r     <- raster::stack(path)
 cells <- raster::cellFromXY(r, xy)
 data  <- raster::extract(r, cells)

，x/y 坐标的矩阵在哪里，xy您可以使用以下内容从 csv 中轻松获得：

 as.matrix(csv_read(cities.path))

高温高压

于 2017-06-16T20:02:46.683 回答