Как извлечь ячейки сетки модели из файла NetCDF, соответствующего определенным широтам и долготам, и показать оба на одном графике?

Question

Я работаю с R только последние три недели, поэтому я очень новичок. В настоящее время я работаю с климатическими данными NetCDF по всей Новой Англии (путь). У меня также есть файл.csv с координатами конкретных городов, которые мы хотим посмотреть (towns.path). Мне удалось извлечь временной годовой временной ряд и тренд из ячеек сетки модели, соответствующих конкретным городам, из файла.csv. Проблема была в том, что я смог нанести на карту важные городские станции из моего файла.csv со средними годовыми данными.

Когда я запускаю строку сценария 'ave_annual_cities <- extract (year_ave_stack, towns.points, df = T)', я получаю график широты / долготы с точками моих важных городов. Когда я запускаю "plot (координировать (towns.points))" в консоли, я снова получаю график широты / долготы с моими важными точками городов, однако он стоит отдельно, отдельно от моего графика ave_annual_cities. Когда я запускаю 'levelplot (subset (year_ave_stack, 1), margin = F) + layer (sp.points (towns.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)

Надеюсь, это имеет смысл.

Answer 1

Если ваши входные данные правильно считываются как многоканальный растр, вы можете прочитать netcdf как растр, используя:

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

, где xy это матрица координат x/y, которую вы можете легко получить из вашего CSV, используя что-то в строках:

 as.matrix(csv_read(cities.path))

НТН