transformation .nc en raster

Postez ici vos questions, réponses, commentaires ou suggestions - Les sujets seront ultérieurement répartis dans les archives par les modérateurs

Modérateur : Groupe des modérateurs

Laëtitia VIBERT
Messages : 42
Enregistré le : 16 Avr 2019, 10:13

transformation .nc en raster

Messagepar Laëtitia VIBERT » 24 Juin 2019, 07:06

Bonjour,
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}
}
L.

Retourner vers « Questions en cours »

Qui est en ligne

Utilisateurs parcourant ce forum : Aucun utilisateur enregistré et 1 invité