気候変数 (気温と降水量) の相関マップを生成したい。私が使用しているパッケージは、ncdf の更新版である「ncdf4」です。私が使用しているスクリプトはncdfによるものですが、ncdf4に従ってコードを実行することはできましたが、実行中にいくつかの問題に直面しています
pre_bl<-pre[west:east, south:north,] # which says subscript out of bounds any solution for this line. I have attached my code herewith this post.
options(jupyter.plot_mimetypes = 'image/png')
setwd("D:/giovanni/correlation/dir")
# Install required Packages
library(chron);
library(RColorBrewer);
library(lattice);
library(ncdf4);
library(maptools);
data(wrld_simpl);
library(fields);
# Open data
# Both datasets are at monthly resolution
pre.nc <- nc_open("near-surface temperature maximum (degrees Celsius).nc")
sst.nc <- nc_open("precipitation(mmpermonth).nc")
# see variables and dimension names
pre.nc
sst.nc
pre <- ncvar_get(pre.nc, "tmx")
sst<-ncvar_get(sst.nc, "pre")
# Next, get latitude and longitude of both variables.
sst.nc$dim$lon$vals -> lon_sst
sst.nc$dim$lat$vals -> lat_sst
pre.nc$dim$lon$vals -> lon_pre
pre.nc$dim$lat$vals -> lat_pre
# Let's close sst.nc and pre.nc
nc_close(sst.nc)
nc_close(pre.nc)
# We can have a look at the latitude and longitude values.
lat_pre[1:80]
lon_pre[1:80]
west <- lon_pre[min(which(lon_pre>60.25))] #
east <- lon_pre[max(which(lon_pre<99.75))] #
south <- lat_pre[min(which(lat_pre>0.25))] #
north <- lat_pre[max(which(lat_pre<39.75))] #
#code worked upto this line
pre_bl<-pre[west:east, south:north,]
dim(pre_bl)
# Now, let's area average it.
pre_areaAve<-apply(pre_bl,3,mean,na.rm=TRUE)
# Now, we can calculate monthly average rainfall over the region to see when the main rainy season is.
monthly_ave<-t(matrix(pre_areaAve,12,33))
monthly_ave<-apply(monthly_ave,2,mean)
# Let's plot the monthly average rainfall to see the main rainy season.
labels<-c("Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec",sep = " ")
title='Monthly Average rainfall'
ylab1="mm/month"
plot(monthly_ave,type='o',col="blue",lwd=3,
ylab=ylab1,main=title,
xaxt = "n", xlab =" ")
axis(1, labels = FALSE)
text(1:12, par("usr")[1] +(min(monthly_ave)+max(monthly_ave))/2, srt = 45, adj = 1,
labels = labels, xpd = TRUE)
sep_pre<-matrix(pre_areaAve,12,33)[9,]
dim(sst)
sep_index <-c()
for (i in 1:33){
sep_index<-c(sep_index,9+(i-1)*12)
}
# Now, we can extract sst in September using sep_index.
sep_sst<-sst[ , ,sep_index]
dim(sep_sst) # we are considering the month of September over 33 years
corr_pcp_sst<-matrix(0,72,36) # created a container
for (i in 1:72){
for (j in 1:36){
corr_pcp_sst[i,j]<-cor(sep_sst[i,j,],sep_pre)
}
}
lat_sst # let's see the latitude values to select region for map display.
# Ok, let's generate the correlation map from 42.5 S to 42.5 N.
mypalette<-c("#67000d","#a50f15","#cb181d","#ef3b2c","#fb6a4a","#fc9272",
"#fcbba1","#fee0d2","#fff5f0")
image.plot(lon_sst,lat_sst[10:27],corr_pcp_sst[ ,10:27], col = mypalette,
xlab ="Longitude", ylab="Latitude",
main ="Correlation between Blue Nile rainfall and SST in September")
plot(wrld_simpl,add=TRUE) # overlay world map