je cherche à 'rasteriser' des données de concentrations de chlorophylle-a (fichiers .nc) sur une zone d'étude, et j'obtiens le message d'erreur suivant, au bout d'un très long moment :
Code : Tout sélectionner
error in .local(x,values,...)
values must be numeric, integer or logical
Problème lié aux valeurs NA ?
Je pensais au début que cela venait de mon ordinateur portable qui n'était pas assez puissant, mais après avoir fait planter de la même manière mon PC au bureau, j'ai des doutes.
Je ne trouve pourtant pas d'erreur dans mon script...
Un immense merci à qui pourra m'éclairer !
Données :
chlorophylle-a: https://oceandata.sci.gsfc.nasa.gov/MOD ... or_a/2019/
(j'utilise les fichiers de 2017 à 2019, mais un seul est nécessaire pour tester le script)
script :
Code : Tout sélectionner
###convert daily sea ice concentration ifremer data to raster in polar stereo coordinates
rm(list=ls())
library(ncdf4) #read netcf file
library(raster) #work with GIS / raster data
library(sp) #work with GIS
library(rgdal) #work with GIS
library(tiff)
library(chron)
library(lattice)
library(RColorBrewer)
path_file_save <- "my_directory" #where to save raster at the end
path_file_AQUA <- "my_directory" #where to find data
setwd(path_file_AQUA)
list <- list.files(pattern = "chlor_a")
list2 <- list.files(pattern = ".x.nc")
for (i in 1:length(list)) {
setwd(path_file_AQUA)
ncfname_Chla <- list[i]
nc_Chla <- nc_open(ncfname_Chla,write=TRUE)
# print(nc_Chla)
dname <- "chlor_a"
# file date
ncfname_Chla
date_Chla <- gsub(".L3m_8D_CHL_chlor_a_4km.nc", "", ncfname_Chla)
year1 <- as.numeric(substr(date_Chla, 2, 5))
year2 <- as.numeric(substr(date_Chla, 9, 12))
day1 <- as.numeric(substr(date_Chla, 6, 8))
day2 <- as.numeric(substr(date_Chla, 13, 15))
origin <- as.vector(c(01,01,year1))
day1 <- month.day.year(day1, origin)
day2 <- month.day.year(day2,origin)
if (day1$month < 10) {day1$month <- paste("0",day1$month, sep="")}
if (day1$day < 10) {day1$day <- paste("0",day1$day, sep="")}
if (day2$month < 10) {day2$month <- paste("0",day2$month, sep="")}
if (day2$day < 10) {day2$day <- paste("0",day2$day, sep="")}
date1 = paste(day1[3],day1[1],day1[2], sep = "")
date2 = paste(day2[3],day2[1],day2[2], sep = "")
date_Chla = paste(date1,date2, sep = "_")
# get coordinates
lon <- ncvar_get(nc_Chla,"lon")
nlon <- dim(lon)
lat <- ncvar_get(nc_Chla,"lat")
nlat <- dim(lat)
# head(lon) ;head(lat)
# print(c(nlon,nlat))
#get Chlorophyll-a data
Chla <- ncvar_get(nc_Chla, dname)
fillvalue <- ncatt_get(nc_Chla,dname,"_FillValue")
Chla[Chla == fillvalue$value] <- NA
#create a df with all coordinates and the Chlor_a concentration
lonlat <- as.matrix(expand.grid(lon,lat))
dim(lonlat) #dim(lon)*dim(lat)
Chla_vec <- as.vector(Chla)
length(Chla_vec)
Chla_all <- data.frame(cbind(lonlat,Chla_vec))
names(Chla_all) <- c("lon","lat",dname)
# define study area
#coordinate of study site
LON_ANT = 140
LAT_ANT = -66.7
Chla_ANT <- Chla_all[which(Chla_all$lon >= (LON_ANT-5) & Chla_all$lon <= (LON_ANT+5)
& Chla_all$lat >= (LAT_ANT-5) & Chla_all$lat <= (LAT_ANT+5)),]
# table(is.na(Chla_ANT$chlor_a)) #NA may occur because of sea-ice, and some periods may be completely unavailable
if (!is.na(mean.default(Chla_ANT$chlor_a, na.rm = T))) {
#rasterize the data
#define coordinates wgs84 and then calculate new ones in UTM or polar stereo to get distance in meters & create raster
setwd(path_file_save)
pts=Chla_ANT[,-3]
coordinates(pts) <- ~lon+lat
proj4string(pts)=CRS("+init=epsg:4326") # set it to lat-long
pts = spTransform(pts,CRS("+init=epsg:3031")) #UTM 54S for DDU
rast <- raster(ncol = 240, nrow = 240)
extent(rast) <- extent(pts)
Chla.rast <- rasterize(pts, rast, Chla_ANT$chlor_a) #le problème vient de là, peut-être dû aux values NA??
proj4string(Chla.rast)=CRS("+init=epsg:3031") # set it to lat-long
setwd(path_file_save)
raster_name <- paste(date_Chla,"Chla_con.tif", sep ="_")
writeRaster(Chla.rast, raster_name, format="GTiff",overwrite=T) #know which projection you are
} else {next}
}