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

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

Когда я запускаю строку сценария «ave_annual_cities‹ - extract (Annual_ave_stack, cities.points, df = T) », я получаю график широты / долготы с точками моих важных городов. Когда я запускаю «plot (координаты (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)

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