Save daily data for each year for 4d array

58 Views Asked by At

I have precipitation data and potential evapotranspiration (PET) data for a period of 5 years. I am undertaking a simple calculation of precipitation/PET for each longitude and latitude within an array. I then want to save this daily data for each month of each year to a netcdf file so I can analyse it in a GIS. However, I have tried many different ways but cannot get it to work.

I am trying to use:

lats = 712
lons = 584

# Initialize an empty 4D array for daily petvals
  petvals_daily = array(0., c(lons, lats,5,365)) #lons, lats, years, days

#Load climate data
  np=nc_open(paste("TAMSAT_", y, "_regrid.nc", sep=''))  

# open precipitation files

  pet=nc_open(paste("cru_ts4.06.pet.dat_remap_", y, "_01_temp_366_final.nc", sep=''))  

# open PET files
  
  P = ncvar_get(np, 'rfe_filled')
  petvalues = ncvar_get(pet,'pet')

  ts=dim(P)[[3]] #set petvals array coordinates the same as precipitation coordinates
  
  # initialize petvals
  petvals = array(0., c(lons, lats, ts))
  
  # read coordinates 
  LON = ncvar_get(np, 'lon')  # set longitude coordinates to same as .nc file with precipitation data
  LAT = ncvar_get(np, 'lat')  # set latitude coordinates to same as .nc file with precipitation data
  ts=dim(P)[[3]] #set petvals array coordinates the same as precipitation coordinates
    
  
  for (i in 1:lats) { 
    for (j in 1:lons) {
      petvals[j, i, ] = petval(start = ss, end = ee, precip = P[j, i, ], PETv = petvals[j, i, ])
    }
  }

# Save daily petvals for each year
  for (d in 1:365) {
    dlDD = format(daterange, '%d') == DD[d]
  petvals_daily[,,y-1982,d] = petvals[,,dlDD]
  }

But this leads to the error:

Error in [<-(*tmp*, , , y - 1982, d, value = petvals[, , dlDD]) : subscript out of bounds

I have tried to modify the subscript placement as well as the structure but it has not worked..

1

There are 1 best solutions below

4
Jonathan Tinsley On

So I answered based on @Onyambu guidance in the comments:

instead of:

# Save daily petvals for each year
  for (d in 1:365) {
    my_years_daily = format(daterange, '%d') == DailyData[d]
  petvals_daily[,,y-1982,d] = petvals[,,my_years_daily]
  }

I needed to simply do:

petvals_daily[,,,] = petvals

Thanks to Onyambu for the guidance!