How to extract pixel values of a polygon shape, using R for annual (365 files per year), daily NetCDF raster files?

You would be much better served by thinking through your analysis and keeping it all in R. It is well know that Excel is very limited in statistical analysis and performing task such as aggregated statistical summaries is very simple, and repeatable, in R. The analytical tools that you have on hand in R, for spatial and statistical analysis, are far superior to attempting to perform an analysis in a spreadsheet type software.

First, I would not recommend "extracting" your time-series data to a polygon per se but, rather keeping it in a raster format masked to the polygon. By extracting to a polygon you introduce unnecessary difficulties into your analysis, and potentially lose the "pixel-level" information. You also make your analysis much less tractable.

Here is a quick run-through example of what you are asking. First, add the required libraries

library(raster)
library(exactextractr)
library(sf)

Now we create some example data that is a stack of 50, 100-row by 100 -column rasters with random values assigned. We also create a vector of dates to work with, that theoretically correspond to our stack of rasters. You would use the raster stack or brick function to read in your raster data into a single object representing the time-series.

i=100;j=100;n=50
r <- do.call(stack, replicate(n,raster(matrix(runif(i*j), i, j))))
  ( dates <- seq(as.Date("2000/1/1"), by = "month", length.out = n) )

Here we create a polygon that would represent an area that we may want to subset. I create an sp class object but you can use the sf st_read function to read in something like a shapefile.

poly <- as(extent(0.1045409,0.9528445,0.1265942,0.8459615), "SpatialPolygons") 
  plot(r[[1]])
    plot(poly,add=TRUE)

Here is where you can simply mask your raster time-series to the polygon AOI. This allows one to perform operations on the time-series in very simple ways.

r.mask <- mask(r, poly) plot(r.mask[[1:12]])

For instance, we can take the pixel means across the time-series as such

r.mean <- calc(r.mask, mean)    
  plot(r.mean)

Or, the mean of a subset of dates (ie., only 2002).

r.mean2002 <- calc(r.mask[[grep("2002", dates)]], mean)

These are very simple examples and for each resulting summary raster you could then easily assign the results to your polygon. If you are resolutely on the path of "extracting" the raster values to a polygon and exporting the data, I would still suggest masking the data and then leveraging one of the sp classes "SpatialPixelsDataFrame" because you can get pixel center coordinates as well as easily write out results from the @data slot.

Here we coerce our raster stack to an sp class, assign the date names, and add columns for the [x,y] coordinates. You can then use write.csv to export the results (a data.frame in the @data slot) to a comma-separated flat file. The columns represent each raster (date) and the rows are the associated pixel values.

r.sp <- as(r, "SpatialPixelsDataFrame")
  names(r.sp) <- as.character(dates)
    r.sp@data <- data.frame(coordinates(r.sp), r.sp@data) 
write.csv(r.sp@data, "my_results.csv") 

I do feel compelled to demonstrate a method for "extracting" the data to a polygon but, in your case, it adds an additional processing step. This would be much more relevant if you had several polygons. You can use the exact_extract function, from the exactextractr package, to actually extract the raster values for the polygon. Please note that there is no efficient way to include [x,y] pixel coordinates so, you just get the pixel values.

r.values <- exact_extract(r, as(poly, "sf"))[[1]]
  names(r.values)[-ncol(r.values)] <- as.character(dates)
    str(r.values)