diff --git a/DESCRIPTION b/DESCRIPTION index e0bb552..8a42251 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,26 +1,26 @@ Package: NASAaccess -Version: 3.4.3 -Date: 2023-04-19 +Version: 4.0.0 +Date: 2024-01-11 Type: Package Title: Downloading and Reformatting Tool for NASA Earth Observation Data Products Authors@R: c(person("Ibrahim", "Mohammed", role = c("aut", "cre"), comment = c(ORCID = "0000-0002-6542-319X"), email = "ibrahim.mohammed@nasa.gov")) Maintainer: Ibrahim Mohammed Depends: R (>= 4.0) -Description: The package assumes that users have already set up a registration account(s) with Earthdata login as well as authorizing NASA GES DISC data access. Please refer to for further details. Creating the '.netrc' file at the user machine Home directory and storing the user NASA GES DISC logging information in it is done automatically to execute the NASAaccess package commands. The GES DISC user registration access logging information will be processed by masking in the terminal on any major OS. Without providing GES DISC user registration access logging information, the user will receive 'You need to provide your Earthdata GES DISC login to proceed…' message. +Description: Users need to set up a registration account(s) with Earthdata login as well as authorizing NASA GES DISC data access to generate gridded ascii tables of climate and earth observation remote sensing data. Please refer to for further details. Creating the '.netrc' file at the user machine Home directory and storing the user NASA GES DISC logging information in it is done automatically to execute the NASAaccess package commands. The GES DISC user registration access logging information will be processed by masking in the terminal on any major OS. Without providing GES DISC user registration access logging information, the user will receive 'You need to provide your Earthdata GES DISC login to proceed…' message. URL: https://github.com/nasa/NASAaccess BugReports: https://github.com/nasa/NASAaccess/issues -Imports: ncdf4, stats, shapefiles, rgeos, maptools, httr, stringr, rgdal, XML, sp, utils, raster, methods, getPass +Imports: ncdf4, stats, httr, stringr, XML, utils, methods, getPass, terra License: file LICENSE Encoding: UTF-8 Language: en-US -RoxygenNote: 7.2.3 +RoxygenNote: 7.3.0 Suggests: markdown, rmarkdown, knitr, - ggmap, - sf, ggplot2, - codetools + codetools, + leaflet, + tidyterra VignetteBuilder: knitr diff --git a/NAMESPACE b/NAMESPACE index 274171f..cd2026a 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -10,19 +10,9 @@ export(NEX_GDDP_CMIP6) import(XML) import(getPass) import(httr) -import(maptools) import(methods) import(ncdf4) -import(rgdal) -import(rgeos) -import(shapefiles) -import(sp) import(stringr) +import(terra, except = c('tail','head')) import(utils) -importFrom(raster,cellFromPolygon) -importFrom(raster,extract) -importFrom(raster,raster) -importFrom(raster,rowColFromCell) -importFrom(raster,xyFromCell) -importFrom(rgeos,gCentroid) importFrom(stats,na.exclude) diff --git a/NEWS.md b/NEWS.md index ddec979..77ef55f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,10 @@ +# NASAaccess 4.0.0 +* Removing raster, rgeos, regal, maptools and transition to terra! +* Fixing CMIP6 server access + +# NASAaccess 3.4.3 +* Adding DOI (10.5281/zenodo.8422392). + # NASAaccess 3.4.3 * Description edits. diff --git a/R/GLDASpolyCentroid.R b/R/GLDASpolyCentroid.R index 1f818b2..d2aaa3c 100644 --- a/R/GLDASpolyCentroid.R +++ b/R/GLDASpolyCentroid.R @@ -1,48 +1,48 @@ -###3/11/19 -#' Generate air temperature input files as well as air temperature stations file from NASA GLDAS remote sensing products. +###1/12/24 +#' Air temperature data from NASA GLDAS on centroid #' #' This function downloads remote sensing data of \acronym{GLDAS} from \acronym{NASA} \acronym{GSFC} servers, extracts air temperature data from grids falling within a specified sub-basin(s) watershed shapefile, and assigns a pseudo air temperature gauge located at the centroid of the sub-basin(s) watershed a weighted-average daily minimum and maximum air temperature data. The function generates tables in a format that \acronym{SWAT} or other rainfall-runoff hydrological model requires for minimum and maximum air temperatures data input. The function also generates the air temperature stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that pseudo grids that correspond to the centroids of the watershed sub-basins. #' @param Dir A directory name to store gridded air temperature and air temperature stations files. -#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). -#' @param DEM A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). +#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'. +#' @param DEM A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'. #' @param start Beginning date for gridded air temperature data. #' @param end Ending date for gridded air temperature data. -#' @details A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASpolyCentroid} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. +#' @details A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASpolyCentroid} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. #' #' The goal of the Global Land Data Assimilation System \acronym{GLDAS} is to ingest satellite and ground-based observational data products, using advanced land surface modeling and data assimilation techniques, in order to generate optimal fields of land surface states and fluxes (Rodell et al., 2004). \acronym{GLDAS} data set used in this function is the \acronym{GLDAS} Noah Land Surface Model L4 3 hourly 0.25 x 0.25 degree V2.1. The full suite of \acronym{GLDAS} data sets are available at \url{https://hydro1.gesdisc.eosdis.nasa.gov/dods/}. The \command{GLDASpolyCentroid} finds the minimum and maximum air temperatures for each day at each grid within the study watershed by searching for minima and maxima over the three hours air temperature data values available for each day and grid. #' -#' The \command{GLDASpolyCentroid} function relies on 'curl' tool to transfer data from \acronym{NASA} servers to a user machine, using HTTPS supported protocol. The 'curl' command embedded in this function to fetch \acronym{GLDAS} netcdf daily global files is designed to work seamlessly given that appropriate logging information are stored in the ".netrc" file and the cookies file ".urs_cookies" as explained in registering with the Earth Observing System Data and Information System. It is imperative to say here that a user machine should have 'curl' installed as a prerequisite to run \command{GLDASpolyCentroid}. +#' The \command{GLDASpolyCentroid} function relies on "curl" tool to transfer data from \acronym{NASA} servers to a user machine, using HTTPS supported protocol. The "curl" command embedded in this function to fetch \acronym{GLDAS} netcdf daily global files is designed to work seamlessly given that appropriate logging information are stored in the ".netrc" file and the cookies file ".urs_cookies" as explained in registering with the Earth Observing System Data and Information System. It is imperative to say here that a user machine should have "curl" installed as a prerequisite to run \command{GLDASpolyCentroid}. #' -#' The \acronym{GLDAS} V2.1 simulation started on January 1, 2000 using the conditions from the \acronym{GLDAS} V2.0 simulation. The \acronym{GLDAS} V2.1 simulation was forced with National Oceanic and Atmospheric Administration \acronym{NOAA}, Global Data Assimilation System \acronym{GDAS} atmospheric analysis fields (Derber et al., 1991), the disaggregated Global Precipitation Climatology Project \acronym{GPCP} precipitation fields (Adler et al., 2003), and the Air Force Weather Agency’s AGRicultural METeorological modeling system \acronym{AGRMET} radiation fields which became available for March 1, 2001 onward. +#' The \acronym{GLDAS} V2.1 simulation started on January 1, 2000 using the conditions from the \acronym{GLDAS} V2.0 simulation. The \acronym{GLDAS} V2.1 simulation was forced with National Oceanic and Atmospheric Administration \acronym{NOAA}, Global Data Assimilation System \acronym{GDAS} atmospheric analysis fields (Derber et al., 1991), the disaggregated Global Precipitation Climatology Project \acronym{GPCP} precipitation fields (Adler et al., 2003), and the Air Force Weather Agency's AGRicultural METeorological modeling system \acronym{AGRMET} radiation fields which became available for March 1, 2001 onward. #' @note #' \command{start} should be equal to or greater than 2000-Jan-01. #' @author Ibrahim Mohammed, \email{ibrahim.mohammed@@nasa.gov} #' @keywords NASA GLDAS Air Temperature #' @return A table that includes points ID, Point file name, Lat, Long, and Elevation information formatted to be read with \acronym{SWAT} or other hydrological model, and #' a scalar of maximum and minimum air temperature gridded data values at a pseudo air temperature grid located at the centroid of each sub-basin within the study watershed provided in ascii format needed by \acronym{SWAT} model or other hydrological model weather inputs. All air temperature tables will be stored at \code{Dir}. -#' @references Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, \cr -#' doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. -#' @references Adler, R. F., G. J. Huffman, A. Chang, R. Ferraro, P.-P. Xie, J. Janowiak, B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. Arkin, and E. Nelkin (2003), The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979–Present), J. Hydrometeorol., 4, 1147-1167, doi:10.1175/1525-7541(2003)004<1147:tvgpcp>2.0.co;2. +#' @references Adler, R. F., G. J. Huffman, A. Chang, R. Ferraro, P.-P. Xie, J. Janowiak, B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. Arkin, and +#' E. Nelkin (2003), The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979–Present), J. Hydrometeorol., 4, 1147-1167, doi:10.1175/1525-7541(2003)004<1147:tvgpcp>2.0.co;2. +#' @references Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. #' @references Rodell, M., P. R. Houser, U. Jambor, J. Gottschalck, K. Mitchell, C.-J. Meng, K. Arsenault, B. Cosgrove, J. Radakovich, M. Bosilovich, J. K. Entin*, J. P. Walker, D. Lohmann, and D. Toll (2004), The Global Land Data Assimilation System, B. Am. Meteorol. Soc., 85, 381-394, doi:10.1175/bams-85-3-381. #' @examples #' #Lower Mekong basin example #' \dontrun{GLDASpolyCentroid(Dir = "./SWAT_INPUT/", watershed = "LowerMekong.shp", #' DEM = "LowerMekong_dem.tif", start = "2015-12-1", end = "2015-12-3")} -#' @import ncdf4 httr stringr rgdal XML utils sp rgeos getPass +#' @import ncdf4 httr stringr utils XML getPass #' @importFrom stats na.exclude -#' @importFrom raster raster cellFromPolygon xyFromCell rowColFromCell extract -#' @importFrom rgeos gCentroid +#' @rawNamespace import(terra, except = c('tail','head')) #' @export + GLDASpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'LowerMekong_dem.tif', start = '2015-12-1', end = '2015-12-3') { - if(file.exists('~/.netrc')==FALSE||file.exists('~/_netrc')==FALSE) + if(file.exists('~/.netrc')==FALSE) { source(system.file("scripts", "netrc.R", package = "NASAaccess")) } - if(file.exists('~/.netrc')==TRUE||file.exists('~/_netrc')==TRUE) + if(file.exists('~/.netrc')==TRUE) { if(length(grep("urs.earthdata.nasa.gov", readLines('~/.netrc')))==!0||length(grep("urs.earthdata.nasa.gov", readLines('~/_netrc')))==!0) { @@ -55,10 +55,11 @@ GLDASpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DE time_period <- seq.Date(from = as.Date(start), to = as.Date(end), by = 'day') # Reading cell elevation data (DEM should be in geographic projection) - watershed.elevation <- raster::raster(DEM) + + watershed.elevation <- terra::rast(DEM) # Reading the study Watershed shapefile - suppressWarnings(polys <- rgdal::readOGR(dsn=watershed,verbose = F)) + polys <- terra::vect(watershed) # SWAT climate 'air temperature' master file name filenametableKEY<-paste(Dir,'temp_Master.txt',sep='') @@ -70,8 +71,9 @@ GLDASpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DE subbasinCentroids <- list() subbasinCentroidsElevation <- list() - subbasinCentroids<-rgeos::gCentroid(polys,byid = TRUE)@coords - subbasinCentroidsElevation<-raster::extract(x=watershed.elevation,y=subbasinCentroids,method='simple') + + suppressWarnings(subbasinCentroids<-terra::crds(terra::centroids(polys))) + subbasinCentroidsElevation<-terra::extract(watershed.elevation, subbasinCentroids) cell.longlatElev<-data.frame(subbasinCentroids,Elev=subbasinCentroidsElevation) names(cell.longlatElev)<-c('LONG','LAT','Elev') @@ -79,17 +81,17 @@ GLDASpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DE #### Begin writing SWAT climate input tables #### Get the SWAT file names and then put the first record date - for(jj in 1:dim(polys@data)[1]) + for(jj in 1:dim(polys)[1]) { if(dir.exists(Dir)==FALSE){dir.create(Dir,recursive = TRUE)} - filenameSWAT[[jj]]<-paste('temp',as.character(polys@data[jj,1]),sep='') + filenameSWAT[[jj]]<-paste('temp',as.character(values(polys[jj,1])),sep='') filenameSWAT_TXT[[jj]]<-paste(Dir,filenameSWAT[[jj]],'.txt',sep='') write(x=format(time_period[1],'%Y%m%d'),file=filenameSWAT_TXT[[jj]]) } #### Write out the SWAT grid information master table - OutSWAT<-data.frame(ID=as.character(polys@data[,1]),NAME=unlist(filenameSWAT),LAT=cell.longlatElev$LAT,LONG=cell.longlatElev$LONG,ELEVATION=cell.longlatElev$Elev) + OutSWAT<-data.frame(ID=as.numeric(unlist(values(polys[,1]))),NAME=unlist(filenameSWAT),LAT=cell.longlatElev$LAT,LONG=cell.longlatElev$LONG,ELEVATION=cell.longlatElev$Elev) utils::write.csv(OutSWAT,filenametableKEY,row.names = F,quote = F) #### Start doing the work! @@ -106,7 +108,7 @@ GLDASpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DE filenames <- XML::readHTMLTable(XML::htmlParse(filenames))[[1]]#getting the subdaily files at each daily URL filenames <- unique(stats::na.exclude(stringr::str_extract(as.character(filenames$Name),'GLDAS.+(.nc4)'))) # Extract the ncdf files - for(ll in 1:length(filenames))# Iterating over each subdaily data file + for(ll in 1:length(filenames))# Iterating over each subdaily data files { # Downloading the file @@ -126,43 +128,45 @@ GLDASpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DE nc.long<-ncdf4::ncvar_get(nc,nc$dim[[lon.dim.number]]) ####getting the x values (latitudes in degrees north) nc.lat<-ncdf4::ncvar_get(nc,nc$dim[[lat.dim.number]]) - data<-array(NA,dim=c(length(nc.lat),length(nc.long),length(filenames))) + + # create a stack + GLDAS<-terra::rast(nrows=length(nc.lat), + ncols=length(nc.long), + nlyrs=length(filenames), + xmin=nc.long[1], + xmax=nc.long[NROW(nc.long)], + ymin=nc.lat[1], + ymax=nc.lat[NROW(nc.lat)], + crs='+proj=longlat +datum=WGS84') } - data[,,ll]<-matrix(as.vector(ncdf4::ncvar_get(nc,nc$var[[33]])),nrow=length(nc.lat),ncol=length(nc.long),byrow=T) - # Reorder the rows - data[,,ll]<-data[nrow(data):1,,ll] + + GLDAS[[ll]] <- as.vector(ncdf4::ncvar_get(nc,nc$var[[33]])) + GLDAS[[ll]] <- terra::flip(GLDAS[[ll]],direction="v") ncdf4::nc_close(nc) - } - # create a stack - GLDAS<-raster::raster(x=data[,,1],xmn=nc.long[1],xmx=nc.long[NROW(nc.long)],ymn=nc.lat[1],ymx=nc.lat[NROW(nc.lat)],crs=sp::CRS('+proj=longlat +datum=WGS84')) - for(j in 1:length(filenames)) - { - rr<-raster::raster(x=data[,,j],xmn=nc.long[1],xmx=nc.long[NROW(nc.long)],ymn=nc.lat[1],ymx=nc.lat[NROW(nc.lat)],crs=sp::CRS('+proj=longlat +datum=WGS84')) - GLDAS<-raster::stack(x=c(GLDAS,rr),quick=T) + } - GLDAS<-raster::dropLayer(GLDAS,1) + ## save time by cropping the world raster to the study DEM - cropGLDAS<-raster::crop(x=GLDAS,y=watershed.elevation) + cropGLDAS<-terra::crop(x=GLDAS,y=watershed.elevation) ### Obtaining subdaily climate values at centroid grids by averaging GLDAS grids data with weights within each subbasin as explained earlier ###daily records for each point - dailytemp<-suppressWarnings(raster::extract(cropGLDAS, polys, weights=TRUE, fun=mean)) - #dailytemp[is.na(dailytemp)] <- '-99.0' #filling missing data + dailytemp<-suppressWarnings(terra::extract(cropGLDAS, polys, weights=TRUE, fun=mean)) ###obtain minimum daily data over the 3 hrs records - mindailytemp<-apply(dailytemp,1,min) + mindailytemp<-apply(dailytemp[,-1],1,min) mindailytemp<-mindailytemp - 273.16 #convert to degree C mindailytemp[is.na(mindailytemp)] <- '-99.0' #filling missing data ###same for maximum daily - maxdailytemp<-apply(dailytemp,1,max) + maxdailytemp<-apply(dailytemp[,-1],1,max) maxdailytemp<-maxdailytemp - 273.16 #convert to degree C maxdailytemp[is.na(maxdailytemp)] <- '-99.0' #filing missing data ### Looping through the GLDAS points and writing out the daily climate data in SWAT format - for(k in 1:dim(polys@data)[1]) + for(k in 1:dim(polys)[1]) { cell.temp.values[[k]]<-paste(maxdailytemp[k],mindailytemp[k],sep=',') write(x=cell.temp.values[[k]],filenameSWAT_TXT[[k]],append=T,ncolumns = 1) @@ -170,7 +174,7 @@ GLDASpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DE #empty memory and getting ready for next day! cell.temp.values<-list() - rm(data,GLDAS) + rm(GLDAS) unlink(x='./temp', recursive = TRUE) } diff --git a/R/GLDASwat.R b/R/GLDASwat.R index 3b34531..d90ff67 100644 --- a/R/GLDASwat.R +++ b/R/GLDASwat.R @@ -1,13 +1,13 @@ -###3/11/19 -#' Generate SWAT air temperature input files as well as air temperature stations file from NASA GLDAS remote sensing products. +###12/28/23 +#' SWAT air temperature data from NASA GLDAS #' #' This function downloads remote sensing data of \acronym{GLDAS} from \acronym{NASA} \acronym{GSFC} servers, extracts air temperature data from grids within a specified watershed shapefile, and then generates tables in a format that \acronym{SWAT} requires for minimum and maximum air temperature data input. The function also generates the air temperature stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that fall within the specified watershed. #' @param Dir A directory name to store gridded air temperature and air temperature stations files. -#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). -#' @param DEM A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). +#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'. +#' @param DEM A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'. #' @param start Beginning date for gridded air temperature data. #' @param end Ending date for gridded air temperature data. -#' @details A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASwat} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. +#' @details A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASwat} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. #' #' The goal of the Global Land Data Assimilation System \acronym{GLDAS} is to ingest satellite and ground-based observational data products, using advanced land surface modeling and data assimilation techniques, in order to generate optimal fields of land surface states and fluxes (Rodell et al., 2004). \acronym{GLDAS} dataset used in this function is the \acronym{GLDAS} Noah Land Surface Model L4 3 hourly 0.25 x 0.25 degree V2.1. The full suite of \acronym{GLDAS} datasets is available at \url{https://hydro1.gesdisc.eosdis.nasa.gov/dods/}. The \command{GLDASwat} finds the minimum and maximum air temperatures for each day at each grid within the study watershed by searching for minima and maxima over the three hours air temperature data values available for each day and grid. #' @@ -20,19 +20,18 @@ #' @keywords NASA GLDAS Air Temperature #' @return A table that includes points ID, Point file name, Lat, Long, and Elevation information formatted to be read with \acronym{SWAT}, and #' a scalar of maximum and minimum air temperature gridded data values at each point within the study watershed in ascii format needed by \acronym{SWAT} model weather inputs will be stored at \code{Dir}. -#' @references Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, \cr -#' doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. #' @references Adler, R. F., G. J. Huffman, A. Chang, R. Ferraro, P.-P. Xie, J. Janowiak, B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. Arkin, and E. Nelkin (2003), The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979–Present), J. Hydrometeorol., 4, 1147-1167, doi:10.1175/1525-7541(2003)004<1147:tvgpcp>2.0.co;2. +#' @references Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. #' @references Rodell, M., P. R. Houser, U. Jambor, J. Gottschalck, K. Mitchell, C.-J. Meng, K. Arsenault, B. Cosgrove, J. Radakovich, M. Bosilovich, J. K. Entin*, J. P. Walker, D. Lohmann, and D. Toll (2004), The Global Land Data Assimilation System, B. Am. Meteorol. Soc., 85, 381-394, doi:10.1175/bams-85-3-381. #' @examples #' #Lower Mekong basin example #' \dontrun{GLDASwat(Dir = "./SWAT_INPUT/", watershed = "LowerMekong.shp", #' DEM = "LowerMekong_dem.tif", start = "2015-12-1", end = "2015-12-3")} -#' @import ncdf4 httr stringr rgdal XML utils sp getPass +#' @import ncdf4 httr stringr utils XML getPass #' @importFrom stats na.exclude -#' @importFrom raster raster cellFromPolygon xyFromCell rowColFromCell extract #' @export + GLDASwat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'LowerMekong_dem.tif', start = '2015-12-1', end = '2015-12-3') { if(file.exists('~/.netrc')==FALSE) @@ -54,10 +53,11 @@ GLDASwat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lowe time_period <- seq.Date(from = as.Date(start), to = as.Date(end), by = 'day') # Reading cell elevation data (DEM should be in geographic projection) - watershed.elevation <- raster::raster(DEM) + watershed.elevation <- terra::rast(DEM) # Reading the study Watershed shapefile - suppressWarnings(polys <- rgdal::readOGR(dsn=watershed,verbose = F)) + + polys <- terra::vect(watershed) # SWAT climate 'air temperature' master file name filenametableKEY<-paste(Dir,'temp_Master.txt',sep='') @@ -95,22 +95,45 @@ GLDASwat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lowe nc.long<-ncdf4::ncvar_get(nc,nc$dim[[lon.dim.number]]) ####getting the x values (latitudes in degrees north) nc.lat<-ncdf4::ncvar_get(nc,nc$dim[[lat.dim.number]]) - #extract data - data<-matrix(as.vector(ncdf4::ncvar_get(nc,nc$var[[33]])),nrow=length(nc.lat),ncol=length(nc.long),byrow=T) - #reorder the rows - data<-data[ nrow(data):1, ] + # create a raster + GLDAS<-terra::rast(nrows=length(nc.lat), + ncols=length(nc.long), + xmin=nc.long[1], + xmax=nc.long[NROW(nc.long)], + ymin=nc.lat[1], + ymax=nc.lat[NROW(nc.lat)], + crs='+proj=longlat +datum=WGS84') + + #fill raster with dummy values + + values(GLDAS) <- 1:ncell(GLDAS) + # Convert raster to points + GLDAS.points <- terra::as.points(GLDAS, na.rm = TRUE) + + # Intersect to keep only points on the shape + GLDAS.points <- GLDAS.points[polys] + + + ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - GLDAS<-raster::raster(x=data,xmn=nc.long[1],xmx=nc.long[NROW(nc.long)],ymn=nc.lat[1],ymx=nc.lat[NROW(nc.lat)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + #obtain cell numbers within the GLDAS raster - cell.no<-raster::cellFromPolygon(GLDAS, polys) + + cell.no <- terra::cells(GLDAS, GLDAS.points)[,2] + #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(GLDAS,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(GLDAS,unlist(cell.no)) - points_elevation<-raster::extract(x=watershed.elevation,y=cell.longlat,method='simple') - study_area_records<-data.frame(ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation) - #sp::coordinates (study_area_records)<- ~x+y - rm(data,GLDAS) + + cell.longlat<-terra::xyFromCell(GLDAS,cell.no) + + + cell.rowCol <- terra::rowColFromCell(GLDAS,cell.no) + + + points_elevation<-terra::extract(x=watershed.elevation,y=cell.longlat,method='simple') + + study_area_records<-data.frame(ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation[,]) + + rm(GLDAS) unlink(x='./temp', recursive = TRUE) @@ -160,15 +183,32 @@ GLDASwat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lowe utils::download.file(quiet = T,method='curl',url=paste(myurl,filenames[ll],sep = ''),destfile = paste('./temp/',filenames[ll],sep = ''), mode = 'wb', extra = '-n -c ~/.urs_cookies -b ~/.urs_cookies -L') # Reading the ncdf file nc<-ncdf4::nc_open( paste('./temp/',filenames[ll],sep = '') ) - data<-matrix(as.vector(ncdf4::ncvar_get(nc,nc$var[[33]])),nrow=length(nc.lat),ncol=length(nc.long),byrow=T) + if(ll==1) + { + GLDAS<-terra::rast(nrows=length(nc.lat), + ncols=length(nc.long), + nlyrs=length(filenames), + xmin=nc.long[1], + xmax=nc.long[NROW(nc.long)], + ymin=nc.lat[1], + ymax=nc.lat[NROW(nc.lat)], + crs='+proj=longlat +datum=WGS84') + } + ### Save the subdaily climate data values in a raster + + GLDAS[[ll]] <- as.vector(ncdf4::ncvar_get(nc,nc$var[[33]])) + + # Reorder the rows - data<-data[ nrow(data):1, ] + GLDAS[[ll]] <- terra::flip(GLDAS[[ll]],direction="v") + + ncdf4::nc_close(nc) - ### Save the subdaily climate data values in a raster - GLDAS<-raster::raster(x=data,xmn=nc.long[1],xmx=nc.long[NROW(nc.long)],ymn=nc.lat[1],ymx=nc.lat[NROW(nc.lat)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + ### Obtaining subdaily climate values at GLDAS grids that has been defined and explained earlier - cell.values<-as.vector(GLDAS)[study_area_records$ID] - SubdailyTemp[[ll]]<-cell.values + + cell.values<-extract(GLDAS[[ll]],study_area_records$ID) + SubdailyTemp[[ll]]<-cell.values[,] } @@ -194,7 +234,7 @@ GLDASwat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lowe #empty memory and getting ready for next day! cell.temp.values<-list() SubdailyTemp<- list() - rm(data,GLDAS) + rm(GLDAS) unlink(x='./temp', recursive = TRUE) } diff --git a/R/GPM_NRT.R b/R/GPM_NRT.R index f20c9cc..5d918b5 100644 --- a/R/GPM_NRT.R +++ b/R/GPM_NRT.R @@ -1,13 +1,13 @@ -###07/05/22 -#' Generate Near Real Time (NRT) rainfall from NASA GPM remote sensing products. +###12/28/23 +#' NASA GPM Near Real Time rainfall data #' #' This function downloads rainfall remote sensing data of \acronym{IMERG} from \acronym{NASA} \acronym{GSFC} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that any hydrological model requires for rainfall data input. The function also generates the rainfall stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that fall within the specified watershed. The minimum latency for this function is one day. #' @param Dir A directory name to store gridded rainfall and rain stations files. -#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). -#' @param DEM A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). +#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'. +#' @param DEM A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'. #' @param start Beginning date for gridded rainfall data. #' @param end Ending date for gridded rainfall data. -#' @details A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDE.06/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products. Units for gridded rainfall data are 'mm'. +#' @details A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDE.06/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products. Units for gridded rainfall data are 'mm'. #' #' \acronym{IMERG} dataset is the GPM Level 3 \acronym{IMERG} *Early* Daily 0.1 x 0.1 deg (GPM_3IMERGDE) derived from the half-hourly \acronym{GPM_3IMERGHHE}. The derived result represents the final estimate of the daily accumulated precipitation. The dataset is produced at the \acronym{NASA} Goddard Earth Sciences (GES) Data and Information Services Center (DISC) by simply summing the valid precipitation retrievals for the day in GPM_3IMERGHHE and giving the result in (mm) \url{https://gpm.nasa.gov/data/directory}. #' @@ -27,9 +27,8 @@ #' #Lower Mekong basin example #' \dontrun{GPM_NRT(Dir = "./INPUT/", watershed = "LowerMekong.shp", #' DEM = "LowerMekong_dem.tif", start = "2022-6-1", end = "2022-6-10")} -#' @import ncdf4 shapefiles rgeos maptools httr stringr rgdal XML utils sp methods getPass +#' @import ncdf4 httr stringr utils XML methods getPass #' @importFrom stats na.exclude -#' @importFrom raster raster cellFromPolygon xyFromCell rowColFromCell extract #' @export @@ -63,10 +62,10 @@ GPM_NRT=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'LowerMekon time_period <- seq.Date(from = as.Date(start), to = as.Date(end), by = 'day') # Reading cell elevation data (DEM should be in geographic projection) - watershed.elevation <- raster::raster(DEM) + watershed.elevation <- terra::rast(DEM) # Reading the study Watershed shapefile - suppressWarnings(polys <- rgdal::readOGR(dsn=watershed,verbose = F)) + polys <- terra::vect(watershed) # weather 'precipitation' master file name filenametableKEY<-paste(Dir,'precipitation','Master.txt',sep='') @@ -106,26 +105,45 @@ GPM_NRT=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'LowerMekon nc.long.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[1]]) ####getting the x values (latitudes in degrees north) nc.lat.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[2]]) - #extract data - data<-ncdf4::ncvar_get(nc,myvarIMERG) - #reorder the rows - data<-data[ nrow(data):1, ] + # create a raster + IMERG<-terra::rast(nrows=length(nc.lat.IMERG), + ncols=length(nc.long.IMERG), + xmin=nc.long.IMERG[1], + xmax=nc.long.IMERG[NROW(nc.long.IMERG)], + ymin=nc.lat.IMERG[1], + ymax=nc.lat.IMERG[NROW(nc.lat.IMERG)], + crs='+proj=longlat +datum=WGS84') + + #fill raster with dummy values + + values(IMERG) <- 1:ncell(IMERG) + ncdf4::nc_close(nc) - ###save the daily weather data values in a raster - IMERG<-raster::raster(x=as.matrix(data),xmn=nc.long.IMERG[1],xmx=nc.long.IMERG[NROW(nc.long.IMERG)],ymn=nc.lat.IMERG[1],ymx=nc.lat.IMERG[NROW(nc.lat.IMERG)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + # Convert raster to points + IMERG.points <- terra::as.points(IMERG, na.rm = TRUE) + + # Intersect to keep only points on the shape + IMERG.points <- IMERG.points[polys] #obtain cell numbers within the IMERG raster - cell.no<-raster::cellFromPolygon(IMERG, polys) + cell.no <- terra::cells(IMERG, IMERG.points)[,2] #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(IMERG,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(IMERG,unlist(cell.no)) - points_elevation<-raster::extract(x=watershed.elevation,y=cell.longlat,method='simple') - study_area_records_IMERG<-data.frame(ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation) - sp::coordinates (study_area_records_IMERG)<- ~x+y - rm(data,IMERG) + cell.longlat<-terra::xyFromCell(IMERG,cell.no) + + + cell.rowCol <- terra::rowColFromCell(IMERG,cell.no) + + + points_elevation<-terra::extract(x=watershed.elevation,y=cell.longlat,method='simple') + + FinalTable<-data.frame(ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation[,]) + + + + rm(IMERG) } - FinalTable = data.frame(sp::coordinates(study_area_records_IMERG),ID=study_area_records_IMERG$ID,row=study_area_records_IMERG$row,col=study_area_records_IMERG$col,Elevation=study_area_records_IMERG$Elevation) + #### Begin writing weather input tables #### Get the file names and then put the first record date @@ -173,15 +191,29 @@ GPM_NRT=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'LowerMekon if(file.exists(paste('./temp/',filenames[ll],sep= ''))==FALSE){utils::download.file(quiet = T,method='curl',url=paste(myurl,filenames[ll],sep = ''),destfile = paste('./temp/',filenames[ll],sep = ''), mode = 'wb', extra = '-n -c ~/.urs_cookies -b ~/.urs_cookies -L')} # Reading the ncdf file nc<-ncdf4::nc_open( paste('./temp/',filenames[ll],sep = '') ) - data<-ncdf4::ncvar_get(nc,myvarIMERG) + if(ll==1) + { + IMERG<-terra::rast(nrows=length(nc.lat.IMERG), + ncols=length(nc.long.IMERG), + xmin=nc.long.IMERG[1], + xmax=nc.long.IMERG[NROW(nc.long.IMERG)], + ymin=nc.lat.IMERG[1], + ymax=nc.lat.IMERG[NROW(nc.lat.IMERG)], + crs='+proj=longlat +datum=WGS84') + } + + ###save the daily weather data values in a raster + values(IMERG) <- ncdf4::ncvar_get(nc,myvarIMERG) + # Reorder the rows - data<-data[ nrow(data):1, ] + + IMERG <- terra::flip(IMERG,direction="v") + ncdf4::nc_close(nc) - ###save the daily weather data values in a raster - IMERG<-raster::raster(x=as.matrix(data),xmn=nc.long.IMERG[1],xmx=nc.long.IMERG[NROW(nc.long.IMERG)],ymn=nc.lat.IMERG[1],ymx=nc.lat.IMERG[NROW(nc.lat.IMERG)],crs=sp::CRS('+proj=longlat +datum=WGS84')) #obtain daily weather values at cells bounded with the study watershed (extract values from a raster) - cell.values<-as.vector(IMERG)[FinalTable$ID] + + cell.values<-extract(IMERG,FinalTable$ID)[,] cell.values[is.na(cell.values)] <- '-99.0' #filling missing data diff --git a/R/GPMployCentroid.R b/R/GPMployCentroid.R index 0f36d5f..1dbbe43 100644 --- a/R/GPMployCentroid.R +++ b/R/GPMployCentroid.R @@ -1,13 +1,13 @@ -###3/11/19 -#' Generate rainfall input files as well as rain station file from NASA GPM remote sensing products. +###8/1/24 +#' NASA GPM rainfall data on centroid #' #' This function downloads rainfall remote sensing data of \acronym{TRMM} and \acronym{IMERG} from \acronym{NASA} \acronym{GSFC} servers, extracts data from grids falling within a specified sub-basin(s) watershed shapefile and assigns a pseudo rainfall gauge located at the centroid of the sub-basin(s) watershed a weighted-average daily rainfall data. The function generates rainfall tables in a format that \acronym{SWAT} or other rainfall-runoff hydrological model requires for rainfall data input. The function also generates the rainfall stations file summary input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those pseudo grids that correspond to the centroids of the watershed sub-basins. #' @param Dir A directory name to store gridded rainfall and rain stations files. -#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). -#' @param DEM A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). +#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'. +#' @param DEM A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'. #' @param start Beginning date for gridded rainfall data. #' @param end Ending date for gridded rainfall data. -#' @details A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. +#' @details A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. #' #' \acronym{IMERG} dataset is the GPM Level 3 \acronym{IMERG} *Final* Daily 0.1 x 0.1 deg (GPM_3IMERGDF) derived from the half-hourly GPM_3IMERGHH. The derived result represents the final estimate of the daily accumulated precipitation. The dataset is produced at the \acronym{NASA} Goddard Earth Sciences (GES) Data and Information Services Center (DISC) by simply summing the valid precipitation retrievals for the day in GPM_3IMERGHH and giving the result in (mm) \url{https://gpm.nasa.gov/data/directory}. #' @@ -27,10 +27,8 @@ #' #Lower Mekong basin example #' \dontrun{GPMpolyCentroid(Dir = "./SWAT_INPUT/", watershed = "LowerMekong.shp", #' DEM = "LowerMekong_dem.tif", start = "2015-12-1", end = "2015-12-3")} -#' @import ncdf4 shapefiles rgeos maptools httr stringr rgdal XML utils sp methods getPass +#' @import ncdf4 httr stringr utils XML methods getPass #' @importFrom stats na.exclude -#' @importFrom raster raster cellFromPolygon xyFromCell rowColFromCell extract -#' @importFrom rgeos gCentroid #' @export @@ -61,13 +59,11 @@ GPMpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM time_period <- seq.Date(from = as.Date(start), to = as.Date(end), by = 'day') # Reading cell elevation data (DEM should be in geographic projection) - watershed.elevation <- raster::raster(DEM) + watershed.elevation <- terra::rast(DEM) # Reading the study Watershed shapefile - suppressWarnings(polys <- rgdal::readOGR(dsn=watershed,verbose = F)) + polys <- terra::vect(watershed) - # Adding a generic column "GRIDCODE" - polys@data$GRIDCODE <- seq(1:dim(polys@data)[1]) # SWAT climate 'precipitation' master file name @@ -80,18 +76,18 @@ GPMpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM subbasinCentroidsElevation <- list() - subbasinCentroids<-rgeos::gCentroid(polys,byid = TRUE)@coords - subbasinCentroidsElevation<-raster::extract(x=watershed.elevation,y=subbasinCentroids,method='simple') + suppressWarnings(subbasinCentroids<-terra::crds(terra::centroids(polys))) + subbasinCentroidsElevation<-terra::extract(watershed.elevation, subbasinCentroids) cell.longlatElev<-data.frame(subbasinCentroids,Elev=subbasinCentroidsElevation) names(cell.longlatElev)<-c('LONG','LAT','Elev') #### Begin writing SWAT climate input tables #### Get the SWAT file names and then put the first record date - for(jj in 1:dim(polys@data)[1]) + for(jj in 1:dim(polys)[1]) { if(dir.exists(Dir)==FALSE){dir.create(Dir,recursive = TRUE)} - filenameSWAT[[jj]]<-paste(myvarTRMM,as.character(polys@data$GRIDCODE[jj]),sep='') + filenameSWAT[[jj]]<-paste(myvarTRMM,as.character(values(polys[jj,1])),sep='') filenameSWAT_TXT[[jj]]<-paste(Dir,filenameSWAT[[jj]],'.txt',sep='') #write the data beginning date once! write(x=format(time_period[1],'%Y%m%d'),file=filenameSWAT_TXT[[jj]]) @@ -99,7 +95,7 @@ GPMpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM #### Write out the SWAT grid information master table - OutSWAT<-data.frame(ID=as.character(polys@data$GRIDCODE),NAME=unlist(filenameSWAT),LAT=cell.longlatElev$LAT,LONG=cell.longlatElev$LONG,ELEVATION=cell.longlatElev$Elev) + OutSWAT<-data.frame(ID=as.numeric(unlist(values(polys[,1]))),NAME=unlist(filenameSWAT),LAT=cell.longlatElev$LAT,LONG=cell.longlatElev$LONG,ELEVATION=cell.longlatElev$Elev) utils::write.csv(OutSWAT,filenametableKEY,row.names = F,quote = F) @@ -136,9 +132,8 @@ GPMpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM utils::download.file(quiet = T,method='curl',url=paste(myurl,filenames[ll],sep = ''),destfile = paste('./temp/',filenames[ll],sep = ''), mode = 'wb', extra = '-n -c ~/.urs_cookies -b ~/.urs_cookies -L') # Reading the ncdf file nc<-ncdf4::nc_open( paste('./temp/',filenames[ll],sep = '') ) - data<-ncdf4::ncvar_get(nc,myvarTRMM) - # Reorder the rows - data<-data[ nrow(data):1, ] + + ###evaluate these values one time! if(ll==1 && kk==1) { @@ -146,19 +141,33 @@ GPMpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM nc.long.TRMM<-ncdf4::ncvar_get(nc,nc$dim[[1]]) ####getting the x values (latitudes in degrees north) nc.lat.TRMM<-ncdf4::ncvar_get(nc,nc$dim[[2]]) + + # create a raster + TRMM<-terra::rast(nrows=length(nc.lat.TRMM), + ncols=length(nc.long.TRMM), + xmin=nc.long.TRMM[1], + xmax=nc.long.TRMM[NROW(nc.long.TRMM)], + ymin=nc.lat.TRMM[1], + ymax=nc.lat.TRMM[NROW(nc.lat.TRMM)], + crs='+proj=longlat +datum=WGS84') } + + + ###save the daily weather data values in a raster + values(TRMM) <- ncdf4::ncvar_get(nc,myvarTRMM) ncdf4::nc_close(nc) - ### Save the daily climate data values in a raster - TRMM<-raster::raster(x=as.matrix(data),xmn=nc.long.TRMM[1],xmx=nc.long.TRMM[NROW(nc.long.TRMM)],ymn=nc.lat.TRMM[1],ymx=nc.lat.TRMM[NROW(nc.lat.TRMM)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + # Reorder the rows + TRMM <- terra::flip(TRMM,direction="v") + ## save time by cropping the world raster to the study DEM - cropTRMM<-raster::crop(x=TRMM,y=watershed.elevation) + cropTRMM<-terra::crop(x=TRMM,y=watershed.elevation) ### Obtaining daily climate values at centroid grids by averaging TRMM grids data with weights within each subbasin as explained earlier - cell.values<-suppressWarnings(raster::extract(cropTRMM, polys, weights=TRUE, fun=mean)) + cell.values<-suppressWarnings(terra::extract(cropTRMM, polys, weights=TRUE, fun=mean)) cell.values[is.na(cell.values)] <- '-99.0' #filling missing data ### Looping through the TRMM points and writing out the daily climate data in SWAT format - for(jj in 1:dim(polys@data)[1]) + for(jj in 1:dim(polys)[1]) { - write(x=cell.values[jj],filenameSWAT_TXT[[jj]],append=T,ncolumns = 1) + write(x=cell.values[jj,-1],filenameSWAT_TXT[[jj]],append=T,ncolumns = 1) } unlink(x='./temp', recursive = TRUE) @@ -188,7 +197,6 @@ GPMpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM utils::download.file(quiet = T,method='curl',url=paste(myurl,filenames[ll],sep = ''),destfile = paste('./temp/',filenames[ll],sep = ''), mode = 'wb', extra = '-n -c ~/.urs_cookies -b ~/.urs_cookies -L') # Reading the ncdf file nc<-ncdf4::nc_open( paste('./temp/',filenames[ll],sep = '') ) - data<-ncdf4::ncvar_get(nc,myvarIMERG) ###evaluate these values one time! if(ll==1) { @@ -196,23 +204,29 @@ GPMpolyCentroid=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM nc.long.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[1]]) ####getting the x values (latitudes in degrees north) nc.lat.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[2]]) + IMERG<-terra::rast(nrows=length(nc.lat.IMERG), + ncols=length(nc.long.IMERG), + xmin=nc.long.IMERG[1], + xmax=nc.long.IMERG[NROW(nc.long.IMERG)], + ymin=nc.lat.IMERG[1], + ymax=nc.lat.IMERG[NROW(nc.lat.IMERG)], + crs='+proj=longlat +datum=WGS84') } - # Reorder the rows - data<-data[ nrow(data):1, ] - ncdf4::nc_close(nc) + ###save the daily climate data values in a raster - IMERG<-raster::raster(x=as.matrix(data),xmn=nc.long.IMERG[1],xmx=nc.long.IMERG[NROW(nc.long.IMERG)],ymn=nc.lat.IMERG[1],ymx=nc.lat.IMERG[NROW(nc.lat.IMERG)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + values(IMERG) <- ncdf4::ncvar_get(nc,myvarIMERG) + ncdf4::nc_close(nc) ## save time by cropping the world raster to the study DEM - cropIMERG<-raster::crop(x=IMERG,y=watershed.elevation) + cropIMERG<-terra::crop(x=IMERG,y=watershed.elevation) #Obtaining daily climate values at centroid grids by averaging IMERG grids data with weights within each subbasin as explained earlier - cell.values<-suppressWarnings(raster::extract(cropIMERG, polys, weights=TRUE, fun=mean)) + cell.values<-suppressWarnings(terra::extract(cropIMERG, polys, weights=TRUE, fun=mean)) cell.values[is.na(cell.values)] <- '-99.0' #filling missing data #loop through the grid points to write out the daily climate data in a SWAT format - for(jj in 1:dim(polys@data)[1]) + for(jj in 1:dim(polys)[1]) { - write(x=cell.values[jj],filenameSWAT_TXT[[jj]],append=T,ncolumns = 1) + write(x=cell.values[jj,-1],filenameSWAT_TXT[[jj]],append=T,ncolumns = 1) } } diff --git a/R/GPMswat.R b/R/GPMswat.R index 82f1a7a..41954cb 100644 --- a/R/GPMswat.R +++ b/R/GPMswat.R @@ -1,13 +1,13 @@ -###3/11/19 -#' Generate SWAT rainfall input files as well as rain stations file from NASA GPM remote sensing products. +###8/1/24 +#' SWAT rainfall data from NASA GPM #' #' This function downloads rainfall remote sensing data of \acronym{TRMM} and \acronym{IMERG} from \acronym{NASA} \acronym{GSFC} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that \acronym{SWAT} requires for rainfall data input. The function also generates the rainfall stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that fall within the specified watershed. #' @param Dir A directory name to store gridded rainfall and rain stations files. -#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). -#' @param DEM A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84'). +#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'. +#' @param DEM A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'. #' @param start Beginning date for gridded rainfall data. #' @param end Ending date for gridded rainfall data. -#' @details A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. +#' @details A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. #' #' \acronym{IMERG} dataset is the GPM Level 3 \acronym{IMERG} *Final* Daily 0.1 x 0.1 deg (GPM_3IMERGDF) derived from the half-hourly GPM_3IMERGHH. The derived result represents the final estimate of the daily accumulated precipitation. The dataset is produced at the \acronym{NASA} Goddard Earth Sciences (GES) Data and Information Services Center (DISC) by simply summing the valid precipitation retrievals for the day in GPM_3IMERGHH and giving the result in (mm) \url{https://gpm.nasa.gov/data/directory}. #' @@ -27,9 +27,8 @@ #' #Lower Mekong basin example #' \dontrun{GPMswat(Dir = "./SWAT_INPUT/", watershed = "LowerMekong.shp", #' DEM = "LowerMekong_dem.tif", start = "2015-12-1", end = "2015-12-3")} -#' @import ncdf4 shapefiles rgeos maptools httr stringr rgdal XML utils sp methods getPass +#' @import ncdf4 httr stringr utils XML methods getPass #' @importFrom stats na.exclude -#' @importFrom raster raster cellFromPolygon xyFromCell rowColFromCell extract #' @export @@ -60,10 +59,10 @@ GPMswat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lower time_period <- seq.Date(from = as.Date(start), to = as.Date(end), by = 'day') # Reading cell elevation data (DEM should be in geographic projection) - watershed.elevation <- raster::raster(DEM) + watershed.elevation <- terra::rast(DEM) # Reading the study Watershed shapefile - suppressWarnings(polys <- rgdal::readOGR(dsn=watershed,verbose = F)) + polys <- terra::vect(watershed) # SWAT climate 'precipitation' master file name filenametableKEY<-paste(Dir,myvarTRMM,'Master.txt',sep='') @@ -103,27 +102,45 @@ GPMswat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lower nc.long.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[1]]) ####getting the x values (latitudes in degrees north) nc.lat.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[2]]) - #extract data - data<-ncdf4::ncvar_get(nc,myvarIMERG) - #reorder the rows - data<-data[ nrow(data):1, ] + # create a raster + IMERG<-terra::rast(nrows=length(nc.lat.IMERG), + ncols=length(nc.long.IMERG), + xmin=nc.long.IMERG[1], + xmax=nc.long.IMERG[NROW(nc.long.IMERG)], + ymin=nc.lat.IMERG[1], + ymax=nc.lat.IMERG[NROW(nc.lat.IMERG)], + crs='+proj=longlat +datum=WGS84') + #fill raster with dummy values + + values(IMERG) <- 1:ncell(IMERG) + ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - IMERG<-raster::raster(x=as.matrix(data),xmn=nc.long.IMERG[1],xmx=nc.long.IMERG[NROW(nc.long.IMERG)],ymn=nc.lat.IMERG[1],ymx=nc.lat.IMERG[NROW(nc.lat.IMERG)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + # Convert raster to points + IMERG.points <- terra::as.points(IMERG, na.rm = TRUE) + + # Intersect to keep only points on the shape + IMERG.points <- IMERG.points[polys] + + #obtain cell numbers within the IMERG raster - cell.no<-raster::cellFromPolygon(IMERG, polys) + cell.no <- terra::cells(IMERG, IMERG.points)[,2] #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(IMERG,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(IMERG,unlist(cell.no)) - points_elevation<-raster::extract(x=watershed.elevation,y=cell.longlat,method='simple') - study_area_records_IMERG<-data.frame(ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation) - sp::coordinates (study_area_records_IMERG)<- ~x+y - rm(data,IMERG) + cell.longlat<-terra::xyFromCell(IMERG,cell.no) + + cell.rowCol <- terra::rowColFromCell(IMERG,cell.no) + + points_elevation<-terra::extract(x=watershed.elevation,y=cell.longlat,method='simple') + + #FinalTable.IMERG<-data.frame(IMERG_ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation[,]) + FinalTable.IMERG<-data.frame(IMERG_ID=unlist(cell.no),cell.rowCol,Elevation=points_elevation[,]) + FinalTable.IMERG.vect<-vect(crs='+proj=longlat +datum=WGS84',cell.longlat, + atts=FinalTable.IMERG) + rm(IMERG) unlink(x='./temp', recursive = TRUE) } # The TRMM data grid information - # Use the same dummy date defined above since TRMM has data up to present with less accurancy. The recomendation is to use IMERG data from 2014-03-12 and onward! + # Use the same dummy date defined above since TRMM has data up to present with less accuracy. The recommendation is to use IMERG data from 2014-03-12 and onward! # update my url with TRMM information myurl = paste(paste(url.TRMM.input,year,mon,sep = '/'),'/',sep = '') @@ -146,42 +163,61 @@ GPMswat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lower nc.long.TRMM<-ncdf4::ncvar_get(nc,nc$dim[[1]]) ####getting the x values (latitudes in degrees north) nc.lat.TRMM<-ncdf4::ncvar_get(nc,nc$dim[[2]]) + # create a raster + TRMM<-terra::rast(nrows=length(nc.lat.TRMM), + ncols=length(nc.long.TRMM), + xmin=nc.long.TRMM[1], + xmax=nc.long.TRMM[NROW(nc.long.TRMM)], + ymin=nc.lat.TRMM[1], + ymax=nc.lat.TRMM[NROW(nc.lat.TRMM)], + crs='+proj=longlat +datum=WGS84') + #fill raster with dummy values + + values(TRMM) <- 1:ncell(TRMM) #gettig the climate data data<-ncdf4::ncvar_get(nc,myvarTRMM) - #reorder the rows - data<-data[ nrow(data):1, ] + ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - TRMM<-raster::raster(x=as.matrix(data),xmn=nc.long.TRMM[1],xmx=nc.long.TRMM[NROW(nc.long.TRMM)],ymn=nc.lat.TRMM[1],ymx=nc.lat.TRMM[NROW(nc.lat.TRMM)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + # Convert raster to points + TRMM.points <- terra::as.points(TRMM, na.rm = TRUE) + + # Intersect to keep only points on the shape + TRMM.points <- TRMM.points[polys] + #obtain cell numbers within the TRMM raster - cell.no<-raster::cellFromPolygon(TRMM, polys) + cell.no <- terra::cells(TRMM, TRMM.points)[,2] #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(TRMM,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(TRMM,unlist(cell.no)) - cell.values<-as.vector(TRMM)[unlist(cell.no)] - study_area_records_TRMM<-data.frame(TRMM_ID=unlist(cell.no),cell.longlat,cell.rowCol) - sp::coordinates (study_area_records_TRMM)<- ~x+y - rm(data,TRMM) + cell.longlat<-terra::xyFromCell(TRMM,cell.no) + + + cell.rowCol <- terra::rowColFromCell(TRMM,cell.no) + + #FinalTable.TRMM<-data.frame(TRMM_ID=unlist(cell.no),cell.longlat,cell.rowCol) + FinalTable.TRMM<-data.frame(TRMM_ID=unlist(cell.no),cell.rowCol) + FinalTable.TRMM.vect<-vect(crs='+proj=longlat +datum=WGS84',cell.longlat, + atts=FinalTable.TRMM) + #study_area_records_TRMM<-data.frame(TRMM_ID=unlist(cell.no),cell.longlat,cell.rowCol) + #sp::coordinates (study_area_records_TRMM)<- ~x+y + rm(TRMM) unlink(x='./temp', recursive = TRUE) } # creating a similarity table that connects IMERG and TRMM grids # calculate euclidean distances to know how to connect TRMM grids with IMERG grids - for (i in 1 : nrow(study_area_records_IMERG)) + for (i in 1 : nrow(FinalTable.IMERG)) { - distVec <- sp::spDistsN1(study_area_records_TRMM,study_area_records_IMERG[i,]) + distVec <- terra::distance(FinalTable.TRMM.vect,FinalTable.IMERG.vect[i,]) minDistVec[[i]] <- min(distVec) closestSiteVec[[i]] <- which.min(distVec) } - PointAssignIDs <- methods::as(study_area_records_TRMM[unlist(closestSiteVec),]$TRMM_ID,'numeric') - PointsAssignCol <- methods::as(study_area_records_TRMM[unlist(closestSiteVec),]$col,'numeric') - PointsAssignRow <- methods::as(study_area_records_TRMM[unlist(closestSiteVec),]$row,'numeric') + PointAssignIDs <- methods::as(FinalTable.TRMM[unlist(closestSiteVec),]$TRMM_ID,'numeric') + PointsAssignCol <- methods::as(FinalTable.TRMM[unlist(closestSiteVec),]$X2,'numeric') + PointsAssignRow <- methods::as(FinalTable.TRMM[unlist(closestSiteVec),]$X1,'numeric') - FinalTable = data.frame(sp::coordinates(study_area_records_IMERG),ID=study_area_records_IMERG$ID,row=study_area_records_IMERG$row,col=study_area_records_IMERG$col,Elevation=study_area_records_IMERG$Elevation, + FinalTable = data.frame(x=crds(FinalTable.IMERG.vect)[,1],y=crds(FinalTable.IMERG.vect)[,2],ID=FinalTable.IMERG$IMERG_ID,row=FinalTable.IMERG$X1,col=FinalTable.IMERG$X2,Elevation=FinalTable.IMERG$Elevation, CloseTRMMIndex=PointAssignIDs,Distance=unlist(minDistVec),TRMMCol=PointsAssignCol,TRMMRow=PointsAssignRow) - #### Begin writing SWAT climate input tables #### Get the SWAT file names and then put the first record date for(jj in 1:dim(FinalTable)[1]) @@ -234,14 +270,27 @@ GPMswat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lower utils::download.file(quiet = T,method='curl',url=paste(myurl,filenames[ll],sep = ''),destfile = paste('./temp/',filenames[ll],sep = ''), mode = 'wb', extra = '-n -c ~/.urs_cookies -b ~/.urs_cookies -L') # Reading the ncdf file nc<-ncdf4::nc_open( paste('./temp/',filenames[ll],sep = '') ) - data<-ncdf4::ncvar_get(nc,myvarTRMM) + if(ll==1) + { + TRMM<-terra::rast(nrows=length(nc.lat.TRMM), + ncols=length(nc.long.TRMM), + xmin=nc.long.TRMM[1], + xmax=nc.long.TRMM[NROW(nc.long.TRMM)], + ymin=nc.lat.TRMM[1], + ymax=nc.lat.TRMM[NROW(nc.lat.TRMM)], + crs='+proj=longlat +datum=WGS84') + } + ###save the daily weather data values in a raster + values(TRMM) <- ncdf4::ncvar_get(nc,myvarTRMM) # Reorder the rows - data<-data[ nrow(data):1, ] + TRMM <- terra::flip(TRMM,direction="v") + + ncdf4::nc_close(nc) - ### Save the daily climate data values in a raster - TRMM<-raster::raster(x=as.matrix(data),xmn=nc.long.TRMM[1],xmx=nc.long.TRMM[NROW(nc.long.TRMM)],ymn=nc.lat.TRMM[1],ymx=nc.lat.TRMM[NROW(nc.lat.TRMM)],crs=sp::CRS('+proj=longlat +datum=WGS84')) + ### Obtaining daily climate values at TRMM grids near the IMERG grids that has been defined and explained earlier - cell.values<-as.vector(TRMM)[FinalTable$CloseTRMMIndex] + + cell.values<-extract(TRMM,FinalTable$CloseTRMMIndex)[,] cell.values[is.na(cell.values)] <- '-99.0' #filling missing data ### Looping through the TRMM points and writing out the daily climate data in SWAT format for(jj in 1:dim(FinalTable)[1]) @@ -276,19 +325,27 @@ GPMswat=function(Dir='./SWAT_INPUT/', watershed ='LowerMekong.shp', DEM = 'Lower utils::download.file(quiet = T,method='curl',url=paste(myurl,filenames[ll],sep = ''),destfile = paste('./temp/',filenames[ll],sep = ''), mode = 'wb', extra = '-n -c ~/.urs_cookies -b ~/.urs_cookies -L') # Reading the ncdf file nc<-ncdf4::nc_open( paste('./temp/',filenames[ll],sep = '') ) - data<-ncdf4::ncvar_get(nc,myvarIMERG) + if(ll==1) + { + IMERG<-terra::rast(nrows=length(nc.lat.IMERG), + ncols=length(nc.long.IMERG), + xmin=nc.long.IMERG[1], + xmax=nc.long.IMERG[NROW(nc.long.IMERG)], + ymin=nc.lat.IMERG[1], + ymax=nc.lat.IMERG[NROW(nc.lat.IMERG)], + crs='+proj=longlat +datum=WGS84') + } + ###save the daily weather data values in a raster + values(IMERG) <- ncdf4::ncvar_get(nc,myvarIMERG) # Reorder the rows - data<-data[ nrow(data):1, ] + IMERG <- terra::flip(IMERG,direction="v") + ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - IMERG<-raster::raster(x=as.matrix(data),xmn=nc.long.IMERG[1],xmx=nc.long.IMERG[NROW(nc.long.IMERG)],ymn=nc.lat.IMERG[1],ymx=nc.lat.IMERG[NROW(nc.lat.IMERG)],crs=sp::CRS('+proj=longlat +datum=WGS84')) #obtain daily climate values at cells bounded with the study watershed (extract values from a raster) - cell.values<-as.vector(IMERG)[FinalTable$ID] + cell.values<-extract(IMERG,FinalTable$ID)[,] cell.values[is.na(cell.values)] <- '-99.0' #filling missing data - #This will get me the mean precipiation averaged over the watershed (this is not a weighted mean) - #v <- extract(IMERG, polys,weights=T,cellnumbers=T,df=T,normalizeWeights=T,sp=T,fun='mean') #loop through the grid points to write out the daily climate data in a SWAT format for(jj in 1:dim(FinalTable)[1]) diff --git a/R/NEXGDDP.R b/R/NEXGDDP.R index b83a0ba..0de0c19 100644 --- a/R/NEXGDDP.R +++ b/R/NEXGDDP.R @@ -1,17 +1,17 @@ -###7/4/22 -#' Generate rainfall or air temperature as well as climate input stations file from NASA NEX-GDDP remote sensing climate change data products needed to drive various hydrological models. +###1/10/24 +#' CMIP5 climate data from NASA NEX-GDDP #' #' This function downloads climate change data of rainfall and air temperature from \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP} \acronym{GSFC} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that any hydrological model requires for rainfall or air temperature data input. The function also generates the climate stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected climatological grids that fall within the specified watershed. The \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP} dataset is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 5 \acronym{CMIP5} and across two of the four greenhouse gas emissions scenarios known as Representative Concentration Pathways \acronym{RCPs} (rcp45, rcp85). #' @param Dir A directory name to store gridded climate data and stations files. -#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'). -#' @param DEM A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'). +#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection crs ='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'. +#' @param DEM A study watershed digital elevation model raster in a geographic projection crs ='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'. #' @param start Beginning date for gridded climate data. #' @param end Ending date for gridded climate data. #' @param model A climate modeling center and name from the World Climate Research Programme \acronym{WCRP} global climate projections through the Coupled Model Intercomparison Project 5 \acronym{CMIP5} (e.g., \acronym{IPSL-CM5A-MR} which is Institut Pierre-Simon Laplace \acronym{CM5A-MR} model). #' @param type A flux data type. It's value can be \acronym{'pr'} for precipitation or \acronym{'tas'} for air temperature. #' @param slice A scenario from the Representative Concentration Pathways. It's value can be \acronym{'rcp45'} , \acronym{'rcp85'}, or \acronym{'historical'}. #' -#' @details A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server for \acronym{NEX-GDDP} climate change data products at (\url{https://www.nccs.nasa.gov/services/climate-data-services}). The function uses variable name ('pr') for rainfall in \acronym{NEX-GDDP} data products and variable name ('tas') for \acronym{NEX-GDDP} minimum ('tasmin') and maximum ('tasmax') air temperature data products. The \command{NEX-GDDP} function outputs gridded rainfall data in 'mm' and gridded air temperature (maximum and minimum) data in degrees 'C'. +#' @details A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server for \acronym{NEX-GDDP} climate change data products at (\url{https://www.nccs.nasa.gov/services/climate-data-services}). The function uses variable name ('pr') for rainfall in \acronym{NEX-GDDP} data products and variable name ('tas') for \acronym{NEX-GDDP} minimum ('tasmin') and maximum ('tasmax') air temperature data products. The \command{NEX-GDDP} function outputs gridded rainfall data in 'mm' and gridded air temperature (maximum and minimum) data in degrees 'C'. #' #' \acronym{NEX-GDDP} dataset is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 5 \acronym{CMIP5} (Taylor et al. 2012) and across two of the four greenhouse gas emissions scenarios known as Representative Concentration Pathways \acronym{RCPs} (Meinshausen et al. 2011). The \acronym{CMIP5} \acronym{GCM} runs were developed in support of the Fifth Assessment Report of the Intergovernmental Panel on Climate Change \acronym{IPCC AR5}. This dataset includes downscaled projections from the 21 models and scenarios for which daily scenarios were produced and distributed under \acronym{CMIP5}. #' The Bias-Correction Spatial Disaggregation \acronym{BCSD} method used in generating the \acronym{NEX-GDDP} dataset is a statistical downscaling algorithm specifically developed to address the current limitations of the global \acronym{GCM} outputs (Wood et al. 2002; Wood et al. 2004; Maurer et al. 2008; Thrasher et al. 2012). The \acronym{NEX-GDDP} climate projections is downscaled at a spatial resolution of 0.25 degrees x 0.25 degrees (approximately 25 km x 25 km). The \command{NEX_GDDP_CMIP5} downscales the \acronym{NEX-GDDP} data to grid points of 0.1 degrees x 0.1 degrees following nearest point methods described by Mohammed et al. (2018). @@ -39,9 +39,8 @@ #' \dontrun{NEX_GDDP_CMIP5(Dir = "./INPUT/", watershed = "LowerMekong.shp", #' DEM = "LowerMekong_dem.tif", start = "2060-12-1", end = "2060-12-3", #' model = 'IPSL-CM5A-MR', type = 'pr', slice = 'rcp85')} -#' @import ncdf4 shapefiles rgeos maptools httr stringr rgdal XML utils sp methods getPass +#' @import ncdf4 httr stringr utils XML methods getPass #' @importFrom stats na.exclude -#' @importFrom raster raster cellFromPolygon xyFromCell rowColFromCell extract #' @export @@ -72,11 +71,11 @@ NEX_GDDP_CMIP5=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low # Constructing time series based on start and end input days! time_period <- seq.Date(from = as.Date(start), to = as.Date(end), by = 'day') # Reading cell elevation data (DEM should be in geographic projection) - watershed.elevation <- raster::raster(DEM) + watershed.elevation <- terra::rast(DEM) # Reading the study Watershed shapefile - suppressWarnings(polys <- rgdal::readOGR(dsn=watershed,verbose = F)) + polys <- terra::vect(watershed) # To address missing parameters in projection strings - suppressWarnings(polys <- sp::spTransform(polys,CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))) + polys <- terra::project(polys, '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') # Grid climate master file name filenametableKEY<-paste(Dir,type, 'Grid_Master.txt',sep='') # Creating empty lists @@ -112,23 +111,38 @@ NEX_GDDP_CMIP5=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low nc.long.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[1]]) ####getting the x values (latitudes in degrees north) nc.lat.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[2]]) + # create a raster + IMERG<-terra::rast(nrows=length(nc.lat.IMERG), + ncols=length(nc.long.IMERG), + xmin=nc.long.IMERG[1], + xmax=nc.long.IMERG[NROW(nc.long.IMERG)], + ymin=nc.lat.IMERG[1], + ymax=nc.lat.IMERG[NROW(nc.lat.IMERG)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') #extract data - data<-ncdf4::ncvar_get(nc,myvarIMERG) + values(IMERG) <- ncdf4::ncvar_get(nc,myvarIMERG) #reorder the rows - data<-data[ nrow(data):1, ] + IMERG <- terra::flip(IMERG,direction="v") ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - IMERG<-raster::raster(x=as.matrix(data),xmn=nc.long.IMERG[1],xmx=nc.long.IMERG[NROW(nc.long.IMERG)],ymn=nc.lat.IMERG[1],ymx=nc.lat.IMERG[NROW(nc.lat.IMERG)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) + # Convert raster to points + IMERG.points <- terra::as.points(IMERG, na.rm = TRUE) + # Intersect to keep only points on the shape + IMERG.points <- IMERG.points[polys] #obtain cell numbers within the IMERG raster - #cell.no<-raster::cellFromPolygon(IMERG, polys) - suppressWarnings(cell.no<-raster::cellFromPolygon(IMERG, polys)) + cell.no <- terra::cells(IMERG, IMERG.points)[,2] #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(IMERG,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(IMERG,unlist(cell.no)) - points_elevation<-raster::extract(x=watershed.elevation,y=cell.longlat,method='simple') - study_area_records_IMERG<-data.frame(ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation) - sp::coordinates (study_area_records_IMERG)<- ~x+y - rm(data,IMERG) + cell.longlat<-terra::xyFromCell(IMERG,cell.no) + + cell.rowCol <- terra::rowColFromCell(IMERG,cell.no) + + points_elevation<-terra::extract(x=watershed.elevation,y=cell.longlat,method='simple') + + + FinalTable.IMERG<-data.frame(IMERG_ID=unlist(cell.no),cell.rowCol,Elevation=points_elevation[,]) + FinalTable.IMERG.vect<-vect(crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs',cell.longlat, + atts=FinalTable.IMERG) + + rm(IMERG) # The NEX-GDDP data grid information # Use the same dummy date defined above since NEX-GDDP has data from 1950 to 2100. # Using dummy date and file info for a file in the NEX-GDDP dataset @@ -153,46 +167,59 @@ NEX_GDDP_CMIP5=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low nc.lat.NEXGDDP<-ncdf4::ncvar_get(nc,nc$dim[[2]]) #getting the climate data for one data as a dummy matrix catch <- ifelse(isTRUE(type=="pr")==TRUE,type,paste(type,'min',sep="")) - data<-ncdf4::ncvar_get(nc,catch, start = c(1,1,1) , count = c(-1, -1 ,1)) - #transpose the data - data <- raster::t(data) + # create a raster + NEX<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + + ###transpose the data and save the daily climate data values in a raster + values(NEX) <- t(ncdf4::ncvar_get(nc,catch, start = c(1,1,1) , count = c(-1, -1 ,1))) + #reorder the rows - data<-data[ nrow(data):1, ] - ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - NEX<-raster::raster(x=as.matrix(data),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) + NEX <- terra::flip(NEX,direction="v") ###rotate the raster to obtain the longitudes extent -180 to 180 - NEX<-raster::rotate(NEX) + NEX <- terra::rotate(NEX) + ncdf4::nc_close(nc) + # Convert raster to points + NEX.points <- terra::as.points(NEX, na.rm = TRUE) + # Intersect to keep only points on the shape + NEX.points <- NEX.points[polys] #obtain cell numbers within the NEX-GDDP raster - #cell.no<-raster::cellFromPolygon(NEX, polys) - suppressWarnings(cell.no<-raster::cellFromPolygon(NEX, polys)) + cell.no <- terra::cells(NEX, NEX.points)[,2] ##check cell.no to address small watershed if(length(unlist(cell.no))==0) { - #cell.no<-raster::cellFromPolygon(NEX, polys, weights = TRUE)[[1]][,"cell"][1] - suppressWarnings(cell.no<-raster::cellFromPolygon(NEX, polys, weights = TRUE)[[1]][,"cell"][1]) + + + cell.no<-terra::cells(NEX, polys, weights = TRUE)[,2] } #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(NEX,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(NEX,unlist(cell.no)) - study_area_records_NEX<-data.frame(NEX_ID=unlist(cell.no),cell.longlat,cell.rowCol) - sp::coordinates (study_area_records_NEX)<- ~x+y - rm(data,NEX) + cell.longlat<-terra::xyFromCell(NEX,cell.no) + cell.rowCol <- terra::rowColFromCell(NEX,cell.no) + FinalTable.NEX<-data.frame(NEX_ID=unlist(cell.no),cell.rowCol) + FinalTable.NEX.vect<-vect(crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs',cell.longlat, + atts=FinalTable.NEX) + + rm(NEX) unlink(x='./temp', recursive = TRUE) } # creating a similarity table that connects IMERG and NEX-GDDP grids # calculate euclidean distances to know how to connect NEX-GDDP grids with IMERG grids - for (i in 1 : nrow(study_area_records_IMERG)) + for (i in 1 : nrow(FinalTable.IMERG)) { - distVec <- sp::spDistsN1(study_area_records_NEX,study_area_records_IMERG[i,]) + distVec <- terra::distance(FinalTable.NEX.vect,FinalTable.IMERG.vect[i,]) minDistVec[[i]] <- min(distVec) closestSiteVec[[i]] <- which.min(distVec) } - PointAssignIDs <- methods::as(study_area_records_NEX[unlist(closestSiteVec),]$NEX_ID,'numeric') - PointsAssignCol <- methods::as(study_area_records_NEX[unlist(closestSiteVec),]$col,'numeric') - PointsAssignRow <- methods::as(study_area_records_NEX[unlist(closestSiteVec),]$row,'numeric') - FinalTable = data.frame(sp::coordinates(study_area_records_IMERG),ID=study_area_records_IMERG$ID,row=study_area_records_IMERG$row,col=study_area_records_IMERG$col,Elevation=study_area_records_IMERG$Elevation, + PointAssignIDs <- methods::as(FinalTable.NEX[unlist(closestSiteVec),]$NEX_ID,'numeric') + PointsAssignCol <- methods::as(FinalTable.NEX[unlist(closestSiteVec),]$X2,'numeric') + PointsAssignRow <- methods::as(FinalTable.NEX[unlist(closestSiteVec),]$X1,'numeric') + FinalTable = data.frame(x=crds(FinalTable.IMERG.vect)[,1],y=crds(FinalTable.IMERG.vect)[,2],ID=FinalTable.IMERG$IMERG_ID,row=FinalTable.IMERG$X1,col=FinalTable.IMERG$X2,Elevation=FinalTable.IMERG$Elevation, CloseNEXIndex=PointAssignIDs,Distance=unlist(minDistVec),NEXCol=PointsAssignCol,NEXRow=PointsAssignRow) #### Begin writing climate input tables #### Get the climate file names and then put the first record date @@ -227,18 +254,26 @@ NEX_GDDP_CMIP5=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low test3<-file.info(paste('./temp/',filename,sep= ''))$size stopifnot('The NEX GDDP server is temporarily unable to service your request due to maintenance downtime or capacity problems. Please try again later.' = test3 > 6.0e6) nc <- ncdf4::nc_open( paste('./temp/',filename,sep = '') ) - data <- ncdf4::ncvar_get(nc,type, start = c(1,1,dayjuilan) , count = c(-1, -1 ,1)) - #transpose the data - data <- raster::t(data) + # create a raster + NEX<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + ###transpose the data and save the daily climate data values in a raster + values(NEX) <- t(ncdf4::ncvar_get(nc,type, start = c(1,1,dayjuilan) , count = c(-1, -1 ,1))) + #reorder the rows - data<-data[ nrow(data):1, ] + NEX <- terra::flip(NEX,direction="v") + ###rotate the raster to obtain the longitudes extent -180 to 180 + NEX <- terra::rotate(NEX) ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - NEX<-raster::raster(x=as.matrix(data),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) - ###rotate to obtain the longitudes in -180 to 180 - NEX<-raster::rotate(NEX) + ### Obtaining daily climate values at NEX grids near the IMERG grids that has been defined and explained earlier, convert units from kg m^-2 s^-1 to mm day^-1 by multiplying with 86400 (60*60*24) - cell.values<-as.vector(NEX)[FinalTable$CloseNEXIndex]*86400 + + cell.values<-terra::extract(NEX,FinalTable$CloseNEXIndex)[,]*86400 cell.values[is.na(cell.values)] <- '-99.0' #filling missing data } @@ -281,31 +316,47 @@ NEX_GDDP_CMIP5=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low test4<-file.info(paste('./temp/',filename_min,sep= ''))$size stopifnot('The NEX GDDP server is temporarily unable to service your request due to maintenance downtime or capacity problems. Please try again later.' = test4 > 6.0e6) nc_min <- ncdf4::nc_open( paste('./temp/',filename_min,sep = '') ) - data_min <- ncdf4::ncvar_get(nc_min,typemin, start = c(1,1,dayjuilan) , count = c(-1, -1 ,1)) - #transpose the data - data_min <- raster::t(data_min) + # create a raster + NEX_min<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + + ###transpose the data and save the daily climate data values in a raster + values(NEX_min) <- t(ncdf4::ncvar_get(nc_min,typemin, start = c(1,1,dayjuilan) , count = c(-1, -1 ,1))) #reorder the rows - data_min<-data_min[ nrow(data_min):1, ] + NEX_min <- terra::flip(NEX_min,direction="v") + ###rotate the raster to obtain the longitudes extent -180 to 180 + NEX <- terra::rotate(NEX_min) ncdf4::nc_close(nc_min) # Reading the tmax ncdf file test5<-file.info(paste('./temp/',filename_max,sep= ''))$size stopifnot('The NEX GDDP server is temporarily unable to service your request due to maintenance downtime or capacity problems. Please try again later.' = test5 > 6.0e6) nc_max <- ncdf4::nc_open( paste('./temp/',filename_max,sep = '') ) - data_max <- ncdf4::ncvar_get(nc_max,typemax, start = c(1,1,dayjuilan) , count = c(-1, -1 ,1)) - #transpose the data - data_max <- raster::t(data_max) + # create a raster + NEX_max<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + + ###transpose the data and save the daily climate data values in a raster + values(NEX_max) <- t(ncdf4::ncvar_get(nc_max,typemax, start = c(1,1,dayjuilan) , count = c(-1, -1 ,1))) + #reorder the rows - data_max<-data_max[ nrow(data_max):1, ] + NEX_max <- terra::flip(NEX_max,direction="v") + ###rotate the raster to obtain the longitudes extent -180 to 180 + NEX <- terra::rotate(NEX_max) ncdf4::nc_close(nc_max) - ###save the daily climate data values in a raster - NEX_min<-raster::raster(x=as.matrix(data_min),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) - NEX_max<-raster::raster(x=as.matrix(data_max),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) - ###rotate to obtain the longitudes in -180 to 180 - NEX_min<-raster::rotate(NEX_min) - NEX_max<-raster::rotate(NEX_max) + ### Obtaining daily climate values at NEX grids near the IMERG grids that has been defined and explained earlier, convert units to C by substracting 273.16 - cell.values_min<-as.vector(NEX_min)[FinalTable$CloseNEXIndex] - 273.16 #convert to degree C - cell.values_max<-as.vector(NEX_max)[FinalTable$CloseNEXIndex] - 273.16 #convert to degree C + cell.values_min<-terra::extract(NEX_min,FinalTable$CloseNEXIndex)[,] - 273.16 #convert to degree C + cell.values_max<-terra::extract(NEX_max,FinalTable$CloseNEXIndex)[,] - 273.16 #convert to degree C cell.values_min[is.na(cell.values_min)] <- '-99.0' #filling missing data cell.values_max[is.na(cell.values_max)] <- '-99.0' #filling missing data } diff --git a/R/NEXGDDP_CMIP6.R b/R/NEXGDDP_CMIP6.R index 2b8911a..76a8590 100644 --- a/R/NEXGDDP_CMIP6.R +++ b/R/NEXGDDP_CMIP6.R @@ -1,24 +1,24 @@ -###11/11/22 -#' Generate rainfall or air temperature as well as climate input stations file from NASA NEX-GDDP-CMIP6 remote sensing climate change data products needed to drive various hydrological models. +###1/10/24 +#' CMIP6 climate data from NASA NEX-GDDP #' #' This function downloads climate change data of rainfall and air temperature from \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP-CMIP6} \acronym{AMES} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that any hydrological model requires for rainfall or air temperature data input. The function also generates the climate stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected climatological grids that fall within the specified watershed. The \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP-CMIP6} data set is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 6 \acronym{CMIP6} and across two of the four "Tier 1" greenhouse gas emissions scenarios known as Shared Socioeconomic Pathways \acronym{SSPs}. #' @param Dir A directory name to store gridded climate data and stations files. -#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'). -#' @param DEM A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'). +#' @param watershed A study watershed shapefile spatially describing polygon(s) in a geographic projection crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'. +#' @param DEM A study watershed digital elevation model raster in a geographic projection crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'. #' @param start Beginning date for gridded climate data. #' @param end Ending date for gridded climate data. #' @param model A climate modeling center and name from the World Climate Research Programme \acronym{WCRP} global climate projections through the Coupled Model Intercomparison Project 6 \acronym{CMIP6} (e.g., \acronym{MIROC6} which is the sixth version of the Model for Interdisciplinary Research on Climate \acronym{MIROC} model). #' @param type A flux data type. It's value can be \acronym{'pr'} for precipitation or \acronym{'tas'} for air temperature. #' @param slice A scenario from the Shared Socioeconomic Pathways (SSPs). It's value can be \acronym{'ssp126'}, \acronym{'ssp245'}, \acronym{'ssp370'}, \acronym{'ssp585'}, or \acronym{'historical'}. #' -#' @details A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. +#' @details A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. #' The function accesses \acronym{NASA} Goddard Space Flight Center server for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} AMES Research Center server for \acronym{NEX-GDDP-CMIP6} #' climate change data products at \url{https://www.nccs.nasa.gov/services/data-collections/land-based-products/nex-gddp-cmip6}. The function uses variable name ('pr') for rainfall in \acronym{NEX-GDDP-CMIP6} data products and variable name ('tas') for \acronym{NEX-GDDP-CMIP6} minimum ('tasmin') and maximum ('tasmax') #' air temperature data products. The \command{NEX_GDDP_CMIP6} function outputs gridded rainfall data in 'mm' and gridded air temperature (maximum and minimum) data in degrees 'C'. #' #' \acronym{NEX-GDDP-CMIP6} dataset is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 6 \acronym{CMIP6} (Eyring et al. 2016) #' and across the four "Tier 1" greenhouse gas emissions scenarios known as Shared Socioeconomic Pathways \acronym{SSPs} (O'Neil et al. 2016; Meinshausen et al. 2020). The \acronym{CMIP6} \acronym{GCM} runs were developed in support of the Sixth Assessment Report -#' of the Intergovernmental Panel on Climate Change \acronym{IPCC AR6}. This data set includes downscaled projections from the 35 models and scenarios for which daily scenarios were produced and distributed under \acronym{CMIP6}. Please visit \acronym{NCCS} Dataportal - Datashare \url{https://ds.nccs.nasa.gov/thredds2/catalog/AMES/NEX/GDDP-CMIP6/catalog.html} to ensure that the requested model and greenhouse gas emissions scenario (\acronym{SSPs}) is available on the server before you execute your run. +#' of the Intergovernmental Panel on Climate Change \acronym{IPCC AR6}. This data set includes downscaled projections from the 35 models and scenarios for which daily scenarios were produced and distributed under \acronym{CMIP6}. Please visit \acronym{NCCS} Dataportal - Datashare technical note \url{https://www.nccs.nasa.gov/sites/default/files/NEX-GDDP-CMIP6-Tech_Note.pdf} to ensure that the requested model and greenhouse gas emissions scenario (\acronym{SSPs}) is available on the server before you execute your run. #' The Bias-Correction Spatial Disaggregation \acronym{BCSD} method used in generating the \acronym{NEX-GDDP-CMIP6} data set is a statistical downscaling algorithm specifically developed to address the current limitations of the global \acronym{GCM} outputs #' (Wood et al. 2002; Wood et al. 2004; Maurer et al. 2008; Thrasher et al. 2012). The \acronym{NEX-GDDP-CMIP6} climate projections is downscaled at a spatial resolution of 0.25 degrees x 0.25 degrees (approximately 25 km x 25 km). #' The \command{NEX_GDDP_CMIP6} downscales the \acronym{NEX-GDDP-CMIP6} data to grid points of 0.1 degrees x 0.1 degrees following nearest point methods described by Mohammed et al. (2018). @@ -48,9 +48,8 @@ #' \dontrun{NEX_GDDP_CMIP6(Dir = "./INPUT/", watershed = "LowerMekong.shp", #' DEM = "LowerMekong_dem.tif", start = "2060-12-1", end = "2060-12-3", #' model = 'MIROC6', type = 'pr', slice = 'ssp245')} -#' @import ncdf4 shapefiles rgeos maptools httr stringr rgdal XML utils sp methods getPass +#' @import ncdf4 httr stringr utils XML methods getPass #' @importFrom stats na.exclude -#' @importFrom raster raster cellFromPolygon xyFromCell rowColFromCell extract #' @export @@ -68,7 +67,7 @@ NEX_GDDP_CMIP6=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low { url.IMERG.input <- 'https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/' - url.GDDP.input <- 'https://ds.nccs.nasa.gov/thredds/ncss/AMES/NEX/GDDP-CMIP6/' + url.GDDP.input <- 'https://ds.nccs.nasa.gov/thredds/ncss/grid/AMES/NEX/GDDP-CMIP6/' url.catalog.input <- paste('https://ds.nccs.nasa.gov/thredds/catalog/AMES/NEX/GDDP-CMIP6/',model,slice,'catalog.html',sep="/") myvarIMERG <- 'precipitationCal' myvarNAME <- 'climate' @@ -76,13 +75,44 @@ NEX_GDDP_CMIP6=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low if(httr::status_code(GET(url.catalog.input))==200) { #let's get the Variant (e.g., r1i1p1f1) - rr <- readLines(url.catalog.input) + rr <- suppressWarnings(readLines(url.catalog.input)) dump.address <- "/thredds/folder.gif" ; my.pattern <- "[[:alnum:]]{8}" folderlines <- grep(dump.address,rr,value=TRUE) - Variant<- unique(as.vector(stringr::str_match_all(folderlines,my.pattern)[[2]])) + #addressing new website layout + if(length(folderlines)==0) + { + dump.address <- c("catalog.html") + folderlines <- grep(dump.address,rr,value=TRUE) + + } + Variant<- terra::unique(as.vector(stringr::str_match_all(folderlines,my.pattern)[[2]])) rm(dump.address,folderlines,my.pattern,rr) - if(type=='pr'){ftp <- paste(url.GDDP.input,model,'/',slice,'/',Variant,'/',type,'/',type,'_','day','_',model,'_',slice,'_',Variant,'_','gn','_',sep='')} - if(type=='tas'){ftp_min <- paste(url.GDDP.input,model,'/',slice,'/',Variant,'/',type,'min','/',type,'min','_','day','_',model,'_',slice,'_',Variant,'_','gn','_',sep='');ftp_max <- paste(url.GDDP.input,model,'/',slice,'/',Variant,'/',type,'max','/',type,'max','_','day','_',model,'_',slice,'_',Variant,'_','gn','_',sep='')} + + #evaluating id from server (e.g., gr or gn) + if(type == 'pr') + { + typeMOD <- type + + } + + if(type == 'tas') + { + typemin <- paste(type,'min',sep='') + typemax <- paste(type,'max',sep='') + typeMOD <- typemax + } + + fttp<-paste('https://ds.nccs.nasa.gov/thredds/catalog/AMES/NEX/GDDP-CMIP6',model,slice,Variant,typeMOD,'catalog.html',sep='/') + r <- httr::GET(fttp) + filenames <- httr::content(r, "text") + filenames <- XML::readHTMLTable(XML::htmlParse(filenames))[[1]][,1] + my.pattern <- paste(Variant,'.+(.nc)',sep='') + filename.suffix <- unique(stats::na.exclude(stringr::str_extract(as.character(filenames),my.pattern)))[1] + suffix <- str_split(filename.suffix,'_')[[1]][2] + rm(r,filenames,my.pattern,filename.suffix) + + if(type=='pr'){ftp <- paste(url.GDDP.input,model,'/',slice,'/',Variant,'/',type,'/',type,'_','day','_',model,'_',slice,'_',Variant,'_',suffix,'_',sep='')} + if(type=='tas'){ftp_min <- paste(url.GDDP.input,model,'/',slice,'/',Variant,'/',type,'min','/',type,'min','_','day','_',model,'_',slice,'_',Variant,'_',suffix,'_',sep='');ftp_max <- paste(url.GDDP.input,model,'/',slice,'/',Variant,'/',type,'max','/',type,'max','_','day','_',model,'_',slice,'_',Variant,'_',suffix,'_',sep='')} ####Before getting to work on this function do this check on start and end dates if (as.Date(start) >= as.Date('1950-01-01') & as.Date(end) <= as.Date('2100-12-31') & slice == 'ssp126' | slice == 'ssp245' | slice == 'ssp370' | slice == 'ssp585' | slice == 'historical') { @@ -90,11 +120,11 @@ NEX_GDDP_CMIP6=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low # Constructing time series based on start and end input days! time_period <- seq.Date(from = as.Date(start), to = as.Date(end), by = 'day') # Reading cell elevation data (DEM should be in geographic projection) - watershed.elevation <- raster::raster(DEM) + watershed.elevation <- terra::rast(DEM) # Reading the study Watershed shapefile - suppressWarnings(polys <- rgdal::readOGR(dsn=watershed,verbose = F)) + polys <- terra::vect(watershed) # To address missing parameters in projection strings - suppressWarnings(polys <- sp::spTransform(polys,CRSobj = c('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'))) + polys <- terra::project(polys, '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') # Hydrological model climate master file name filenametableKEY<-paste(Dir,type, 'Grid_Master.txt',sep='') # Creating empty lists @@ -130,80 +160,102 @@ NEX_GDDP_CMIP6=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low nc.long.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[1]]) ####getting the x values (latitudes in degrees north) nc.lat.IMERG<-ncdf4::ncvar_get(nc,nc$dim[[2]]) + # create a raster + IMERG<-terra::rast(nrows=length(nc.lat.IMERG), + ncols=length(nc.long.IMERG), + xmin=nc.long.IMERG[1], + xmax=nc.long.IMERG[NROW(nc.long.IMERG)], + ymin=nc.lat.IMERG[1], + ymax=nc.lat.IMERG[NROW(nc.lat.IMERG)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') #extract data - data<-ncdf4::ncvar_get(nc,myvarIMERG) + values(IMERG) <- ncdf4::ncvar_get(nc,myvarIMERG) #reorder the rows - data<-data[ nrow(data):1, ] + IMERG <- terra::flip(IMERG,direction="v") ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - IMERG<-raster::raster(x=as.matrix(data),xmn=nc.long.IMERG[1],xmx=nc.long.IMERG[NROW(nc.long.IMERG)],ymn=nc.lat.IMERG[1],ymx=nc.lat.IMERG[NROW(nc.lat.IMERG)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) + # Convert raster to points + IMERG.points <- terra::as.points(IMERG, na.rm = TRUE) + # Intersect to keep only points on the shape + IMERG.points <- IMERG.points[polys] #obtain cell numbers within the IMERG raster - #cell.no<-raster::cellFromPolygon(IMERG, polys) - suppressWarnings(cell.no<-raster::cellFromPolygon(IMERG, polys)) + cell.no <- terra::cells(IMERG, IMERG.points)[,2] + #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(IMERG,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(IMERG,unlist(cell.no)) - points_elevation<-raster::extract(x=watershed.elevation,y=cell.longlat,method='simple') - study_area_records_IMERG<-data.frame(ID=unlist(cell.no),cell.longlat,cell.rowCol,Elevation=points_elevation) - sp::coordinates (study_area_records_IMERG)<- ~x+y - rm(data,IMERG) + cell.longlat<-terra::xyFromCell(IMERG,cell.no) + cell.rowCol <- terra::rowColFromCell(IMERG,cell.no) + points_elevation<-terra::extract(x=watershed.elevation,y=cell.longlat,method='simple') + FinalTable.IMERG<-data.frame(IMERG_ID=unlist(cell.no),cell.rowCol,Elevation=points_elevation[,]) + FinalTable.IMERG.vect<-vect(crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs',cell.longlat, + atts=FinalTable.IMERG) + rm(IMERG) # The NEX-GDDP-CMIP6 data grid information # Use the same dummy date defined above since NEX-GDDP-CMIP6 has data from 1950 to 2100. # Using dummy date and file info for a file in the NEX-GDDP-CMIP6 dataset # downloading one file if(dir.exists('./temp/')==FALSE){dir.create('./temp/')} - utils::download.file(quiet = T, method = 'curl', url = 'https://ds.nccs.nasa.gov/thredds/ncss/AMES/NEX/GDDP-CMIP6/ACCESS-CM2/ssp585/r1i1p1f1/tasmax/tasmax_day_ACCESS-CM2_ssp585_r1i1p1f1_gn_2015.nc?var=tasmax&disableLLSubset=on&disableProjSubset=on&horizStride=1&time_start=2015-09-01T12%3A00%3A00Z&time_end=2015-09-02T12%3A00%3A00Z&timeStride=1', destfile = paste('./temp/','tasmax_day_ssp585_r1i1p1f1_ACCESS-CM2_2015.nc',sep= ''), mode = 'wb', extra = '-k') - test2<-file.info(paste('./temp/','tasmax_day_ssp585_r1i1p1f1_ACCESS-CM2_2015.nc',sep= ''))$size + utils::download.file(quiet = T, method = 'curl', url = 'https://ds.nccs.nasa.gov/thredds/ncss/grid/AMES/NEX/GDDP-CMIP6/ACCESS-CM2/ssp585/r1i1p1f1/tasmax/tasmax_day_ACCESS-CM2_ssp585_r1i1p1f1_gn_2015.nc?var=tasmax&disableLLSubset=on&disableProjSubset=on&horizStride=1&time_start=2015-09-01T12%3A00%3A00Z&time_end=2015-09-02T12%3A00%3A00Z&timeStride=1', destfile = paste('./temp/','tasmax_day_ssp585_r1i1p1f1_ACCESS-CM2_2015.nc',sep= ''), mode = 'wb', extra = '-k') + test2<-file.info(paste('./temp/','tasmax_day_ssp585_r1i1p1f1_ACCESS-CM2_2015.nc',sep= ''))$size stopifnot('The NEX GDDP server is temporarily unable to service your request due to maintenance downtime or capacity problems. Please try again later.' = test2 > 6.0e6) #reading ncdf file nc<-ncdf4::nc_open( paste('./temp/','tasmax_day_ssp585_r1i1p1f1_ACCESS-CM2_2015.nc',sep = '') ) #since geographic info for all NEX files are the same ###evaluate these values at one time! ###getting the x values (longitudes in degrees east, 0 to +360) so it needed to be converted to -180 to 180) - nc.long.NEXGDDP<-ncdf4::ncvar_get(nc,nc$dim[[3]]) - ####getting the y values (latitudes in degrees north, -90 to +90) - nc.lat.NEXGDDP<-ncdf4::ncvar_get(nc,nc$dim[[2]]) - #getting the climate data - data<-ncdf4::ncvar_get(nc,'tasmax', start = c(1,1,1) , count = c(-1, -1 ,1)) - #transpose the data - data <- raster::t(data) + nc.long.NEXGDDP<-ncdf4::ncvar_get(nc,nc$dim[[1]]) + ####getting the y values (latitudes in degrees north, -60 to +90) + nc.lat.NEXGDDP<-ncdf4::ncvar_get(nc,nc$dim[[3]]) + # create a raster + NEX<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + ###transpose the data and save the daily climate data values in a raster + values(NEX) <- t(ncdf4::ncvar_get(nc,'tasmax', start = c(1,1,1) , count = c(-1, -1 ,1))) + #reorder the rows - data<-data[ nrow(data):1, ] - ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - NEX<-raster::raster(x=as.matrix(data),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) + NEX <- terra::flip(NEX,direction="v") ###rotate the raster to obtain the longitudes extent -180 to 180 - NEX<-raster::rotate(NEX) + NEX <- terra::rotate(NEX) + ncdf4::nc_close(nc) + # Convert raster to points + NEX.points <- terra::as.points(NEX, na.rm = TRUE) + # Intersect to keep only points on the shape + NEX.points <- NEX.points[polys] #obtain cell numbers within the NEX-GDDP raster - #cell.no<-raster::cellFromPolygon(NEX, polys) - suppressWarnings(cell.no<-raster::cellFromPolygon(NEX, polys)) + cell.no <- terra::cells(NEX, NEX.points)[,2] + ##check cell.no to address small watershed if(length(unlist(cell.no))==0) { - #cell.no<-raster::cellFromPolygon(NEX, polys, weights = TRUE)[[1]][,"cell"][1] - suppressWarnings(cell.no<-raster::cellFromPolygon(NEX, polys, weights = TRUE)[[1]][,"cell"][1]) + + cell.no<-terra::cells(NEX, polys, weights = TRUE)[,2] } #obtain lat/long values corresponding to watershed cells - cell.longlat<-raster::xyFromCell(NEX,unlist(cell.no)) - cell.rowCol <- raster::rowColFromCell(NEX,unlist(cell.no)) - study_area_records_NEX<-data.frame(NEX_ID=unlist(cell.no),cell.longlat,cell.rowCol) - sp::coordinates (study_area_records_NEX)<- ~x+y - rm(data,NEX) + cell.longlat<-terra::xyFromCell(NEX,cell.no) + cell.rowCol <- terra::rowColFromCell(NEX,cell.no) + + FinalTable.NEX<-data.frame(NEX_ID=unlist(cell.no),cell.rowCol) + FinalTable.NEX.vect<-vect(crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs',cell.longlat, + atts=FinalTable.NEX) + rm(NEX) } # creating a similarity table that connects IMERG and TRMM grids # calculate euclidean distances to know how to connect TRMM grids with IMERG grids - for (i in 1 : nrow(study_area_records_IMERG)) + for (i in 1 : nrow(FinalTable.IMERG)) { - distVec <- sp::spDistsN1(study_area_records_NEX,study_area_records_IMERG[i,]) + distVec <- terra::distance(FinalTable.NEX.vect,FinalTable.IMERG.vect[i,]) minDistVec[[i]] <- min(distVec) closestSiteVec[[i]] <- which.min(distVec) } - PointAssignIDs <- methods::as(study_area_records_NEX[unlist(closestSiteVec),]$NEX_ID,'numeric') - PointsAssignCol <- methods::as(study_area_records_NEX[unlist(closestSiteVec),]$col,'numeric') - PointsAssignRow <- methods::as(study_area_records_NEX[unlist(closestSiteVec),]$row,'numeric') - FinalTable = data.frame(sp::coordinates(study_area_records_IMERG),ID=study_area_records_IMERG$ID,row=study_area_records_IMERG$row,col=study_area_records_IMERG$col,Elevation=study_area_records_IMERG$Elevation, + PointAssignIDs <- methods::as(FinalTable.NEX[unlist(closestSiteVec),]$NEX_ID,'numeric') + PointsAssignCol <- methods::as(FinalTable.NEX[unlist(closestSiteVec),]$X2,'numeric') + PointsAssignRow <- methods::as(FinalTable.NEX[unlist(closestSiteVec),]$X1,'numeric') + FinalTable = data.frame(x=crds(FinalTable.IMERG.vect)[,1],y=crds(FinalTable.IMERG.vect)[,2],ID=FinalTable.IMERG$IMERG_ID,row=FinalTable.IMERG$X1,col=FinalTable.IMERG$X2,Elevation=FinalTable.IMERG$Elevation, CloseNEXIndex=PointAssignIDs,Distance=unlist(minDistVec),NEXCol=PointsAssignCol,NEXRow=PointsAssignRow) #### Begin writing hydrological model climate input tables #### Get the hydrological model file names and then put the first record date @@ -237,18 +289,25 @@ NEX_GDDP_CMIP6=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low test3<-file.info(paste('./temp/',filename,sep= ''))$size stopifnot('The NEX GDDP server is temporarily unable to service your request due to maintenance downtime or capacity problems. Please try again later.' = test3 > 3.0e6) nc <- ncdf4::nc_open( paste('./temp/',filename,sep = '') ) - data <- ncdf4::ncvar_get(nc,type, start = c(1,1,1) , count = c(-1, -1 ,1)) - #transpose the data - data <- raster::t(data) + # create a raster + NEX<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + + ###transpose the data and save the daily climate data values in a raster + values(NEX) <- t(ncdf4::ncvar_get(nc,type, start = c(1,1,1) , count = c(-1, -1 ,1))) #reorder the rows - data<-data[ nrow(data):1, ] + NEX <- terra::flip(NEX,direction="v") + ###rotate the raster to obtain the longitudes extent -180 to 180 + NEX <- terra::rotate(NEX) ncdf4::nc_close(nc) - ###save the daily climate data values in a raster - NEX<-raster::raster(x=as.matrix(data),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) - ###rotate to obtain the longitudes in -180 to 180 - NEX<-raster::rotate(NEX) + ### Obtaining daily climate values at NEX grids near the IMERG grids that has been defined and explained earlier, convert units from kg m^-2 s^-1 to mm day^-1 by multiplying with 86400 (60*60*24) - cell.values<-as.vector(NEX)[FinalTable$CloseNEXIndex]*86400 + cell.values<-terra::extract(NEX,FinalTable$CloseNEXIndex)[,]*86400 cell.values[is.na(cell.values)] <- '-99.0' #filling missing data ### Looping through the NEX points and writing out the daily climate data in hydrology model format for(jj in 1:dim(FinalTable)[1]) @@ -267,8 +326,7 @@ NEX_GDDP_CMIP6=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low timestart <- time_period[jj] timeend <- timestart + 1 timeyear <- format(timestart,"%Y") - typemin <- paste(type,'min',sep='') - typemax <- paste(type,'max',sep='') + filename_min <- paste(typemin,'_day_',slice,'_',Variant,'_',model,'_',as.character(timestart),'_',as.character(timeend),'.nc',sep = '') filename_max <- paste(typemax,'_day_',slice,'_',Variant,'_',model,'_',as.character(timestart),'_',as.character(timeend),'.nc',sep = '') myurl_min <- paste(ftp_min,timeyear,'.nc?','var=',type,'min','&disableLLSubset=on&disableProjSubset=on&horizStride=1&time_start=',as.character(timestart),'T12%3A00%3A00Z&time_duration=P1D','&timeStride=1',sep = '') @@ -280,31 +338,52 @@ NEX_GDDP_CMIP6=function(Dir='./INPUT/', watershed ='LowerMekong.shp', DEM = 'Low test4<-file.info(paste('./temp/',filename_min,sep= ''))$size stopifnot('The NEX GDDP server is temporarily unable to service your request due to maintenance downtime or capacity problems. Please try again later.' = test4 > 3.0e6) nc_min <- ncdf4::nc_open( paste('./temp/',filename_min,sep = '') ) - data_min <- ncdf4::ncvar_get(nc_min,typemin, start = c(1,1,1) , count = c(-1, -1 ,1)) - #transpose the data - data_min <- raster::t(data_min) + # create a raster + NEX_min<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + + ###transpose the data and save the daily climate data values in a raster + values(NEX_min) <- t(ncdf4::ncvar_get(nc_min,typemin, start = c(1,1,1) , count = c(-1, -1 ,1))) + #reorder the rows - data_min<-data_min[ nrow(data_min):1, ] + NEX_min <- terra::flip(NEX_min,direction="v") + ###rotate the raster to obtain the longitudes extent -180 to 180 + NEX <- terra::rotate(NEX_min) ncdf4::nc_close(nc_min) # Reading the tmax ncdf file test5<-file.info(paste('./temp/',filename_max,sep= ''))$size stopifnot('The NEX GDDP server is temporarily unable to service your request due to maintenance downtime or capacity problems. Please try again later.' = test5 > 3.0e6) nc_max <- ncdf4::nc_open( paste('./temp/',filename_max,sep = '') ) - data_max <- ncdf4::ncvar_get(nc_max,typemax, start = c(1,1,1) , count = c(-1, -1 ,1)) - #transpose the data - data_max <- raster::t(data_max) + # create a raster + NEX_max<-terra::rast(nrows=length(nc.lat.NEXGDDP), + ncols=length(nc.long.NEXGDDP), + xmin=nc.long.NEXGDDP[1], + xmax=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)], + ymin=nc.lat.NEXGDDP[1], + ymax=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)], + crs='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs') + + ###transpose the data and save the daily climate data values in a raster + values(NEX_max) <- t(ncdf4::ncvar_get(nc_max,typemax, start = c(1,1,1) , count = c(-1, -1 ,1))) + + #reorder the rows - data_max<-data_max[ nrow(data_max):1, ] + NEX_max <- terra::flip(NEX_max,direction="v") + ###rotate the raster to obtain the longitudes extent -180 to 180 + NEX <- terra::rotate(NEX_max) ncdf4::nc_close(nc_max) - ###save the daily climate data values in a raster - NEX_min<-raster::raster(x=as.matrix(data_min),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) - NEX_max<-raster::raster(x=as.matrix(data_max),xmn=nc.long.NEXGDDP[1],xmx=nc.long.NEXGDDP[NROW(nc.long.NEXGDDP)],ymn=nc.lat.NEXGDDP[1],ymx=nc.lat.NEXGDDP[NROW(nc.lat.NEXGDDP)],crs=sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs')) - ###rotate to obtain the longitudes in -180 to 180 - NEX_min<-raster::rotate(NEX_min) - NEX_max<-raster::rotate(NEX_max) + ### Obtaining daily climate values at NEX grids near the IMERG grids that has been defined and explained earlier, convert units to C by substracting 273.16 - cell.values_min<-as.vector(NEX_min)[FinalTable$CloseNEXIndex] - 273.16 #convert to degree C - cell.values_max<-as.vector(NEX_max)[FinalTable$CloseNEXIndex] - 273.16 #convert to degree C + + cell.values_min<-terra::extract(NEX_min,FinalTable$CloseNEXIndex)[,] - 273.16 #convert to degree C + cell.values_max<-terra::extract(NEX_max,FinalTable$CloseNEXIndex)[,] - 273.16 #convert to degree C + + cell.values_min[is.na(cell.values_min)] <- '-99.0' #filling missing data cell.values_max[is.na(cell.values_max)] <- '-99.0' #filling missing data diff --git a/README.Rmd b/README.Rmd index 5ffd43f..90f5fc7 100644 --- a/README.Rmd +++ b/README.Rmd @@ -53,7 +53,7 @@ On a local machine the user should have installed the following programs as well - [Installing Rstudio software](https://posit.co/) (OPTIONAL) -- *NASAaccess* R package needs a user registration access with [Earthdata](https://www.earthdata.nasa.gov/). Users should set up a registration account(s) with [Earthdata](https://www.earthdata.nasa.gov/) login as well as well as authorizing [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") [GES DISC](https://disc.gsfc.nasa.gov/) data access. Please refer to for further details. +- *NASAaccess* R package needs a user registration access with [Earthdata](https://www.earthdata.nasa.gov/). Users should set up a registration account(s) with [Earthdata](https://www.earthdata.nasa.gov/) login as well as well as authorizing [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") [GES DISC](https://disc.gsfc.nasa.gov/) data access. Please refer to Data Access section for further details. - [*curl*](https://curl.se/) software . Since Mac users have [*curl*](https://curl.se/) as part of macOS build, Windows users should make sure that their local machines build have [*curl*](https://curl.se/) installed properly. @@ -70,7 +70,7 @@ On a local machine the user should have installed the following programs as well Within the Rstudio help tab the user can verify that the package has been installed and browse the help pages of the various functions of *NASAaccess*. The [GES DISC](https://disc.gsfc.nasa.gov/) user registration access logging information will be processed by masking (i.e., not displaying the lieteral typed text as input) on most but not all platforms. Without providing [GES DISC](https://disc.gsfc.nasa.gov/) user registration access logging information, the user will receive '*You need to provide your Earthdata GES DISC login to proceed...*' message. ### **Conda Package** -- *NASAaccess* conda package needs a user registration access with [Earthdata](https://www.earthdata.nasa.gov/). Users should set up a registration account(s) with [Earthdata](https://www.earthdata.nasa.gov/) login as well as well as authorizing [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") [GES DISC](https://disc.gsfc.nasa.gov/) data access. Please refer to for further details. +- *NASAaccess* conda package needs a user registration access with [Earthdata](https://www.earthdata.nasa.gov/). Users should set up a registration account(s) with [Earthdata](https://www.earthdata.nasa.gov/) login as well as well as authorizing [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") [GES DISC](https://disc.gsfc.nasa.gov/) data access. Please refer to Data Access section for further details. - Creating the *.netrc* file at the user machine *Home* directory and storing the user [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") [GES DISC](https://disc.gsfc.nasa.gov/) logging information in it is done automatically to execute the *NASAaccess* package commands. Accessing data from NASA servers is further explained at [Here](https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+cURL+And+Wget). The [GES DISC](https://disc.gsfc.nasa.gov/) user registration access logging information will be processed by masking in the terminal on any major OS. Without providing [GES DISC](https://disc.gsfc.nasa.gov/) user registration access logging information, the user will receive '*You need to provide your Earthdata GES DISC login to proceed...*' message. @@ -96,3 +96,6 @@ citation(package = 'NASAaccess') ``` +## Reference + +Mohammed, I.N., Bustamante, E.G.R., Bolten, J.D., Nelson, E.J., 2023. Technical note: NASAaccess – a tool for access, reformatting, and visualization of remotely sensed earth observation and climate data. Hydrol. Earth Syst. Sci. 27, 3621-3642, https://doi.org/10.5194/hess-27-3621-2023 diff --git a/README.md b/README.md index 21c3e37..62b0305 100644 --- a/README.md +++ b/README.md @@ -16,8 +16,7 @@ Badge](https://anaconda.org/conda-forge/r-nasaaccess/badges/platforms.svg)](http + [*Ibrahim N. Mohammed*](https://science.gsfc.nasa.gov/sed/bio/ibrahim.mohammed "Ibrahim N. Mohammed") @@ -79,7 +78,8 @@ local machine: authorizing [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") [GES DISC](https://disc.gsfc.nasa.gov/) data access. Please refer to - for further details. + Data Access section + for further details. - [*curl*](https://curl.se/) software . Since Mac users have [*curl*](https://curl.se/) as part of macOS build, Windows users @@ -119,7 +119,8 @@ message. authorizing [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") [GES DISC](https://disc.gsfc.nasa.gov/) data access. Please refer to - for further details. + Data Access section + for further details. - Creating the *.netrc* file at the user machine *Home* directory and storing the user @@ -158,11 +159,10 @@ functionality and capabilities. ``` r citation(package = 'NASAaccess') -#> #> To cite package 'NASAaccess' in publications use: #> -#> Mohammed I (2023). _NASAaccess: Downloading and Reformatting Tool for -#> NASA Earth Observation Data Products_. R package version 3.4.2, +#> Mohammed I (2024). _NASAaccess: Downloading and Reformatting Tool for +#> NASA Earth Observation Data Products_. R package version 4.0.0, #> . #> #> A BibTeX entry for LaTeX users is @@ -170,10 +170,18 @@ citation(package = 'NASAaccess') #> @Manual{, #> title = {{NASAaccess}: Downloading and Reformatting Tool for NASA Earth Observation Data Products}, #> author = {Ibrahim Mohammed}, -#> year = {2023}, +#> year = {2024}, #> institution = {National Aeronautics and Space Administration, Goddard Space Flight Center}, #> address = {Greenbelt, Maryland}, -#> note = {R package version 3.4.2}, +#> note = {R package version 4.0.0}, #> url = {https://github.com/nasa/NASAaccess}, #> } ``` + +## Reference + +Mohammed, I.N., Bustamante, E.G.R., Bolten, J.D., Nelson, E.J., 2023. +Technical note: NASAaccess – a tool for access, reformatting, and +visualization of remotely sensed earth observation and climate data. +Hydrol. Earth Syst. Sci. 27, 3621-3642, + diff --git a/man/GLDASpolyCentroid.Rd b/man/GLDASpolyCentroid.Rd index a7a066a..7ac2733 100644 --- a/man/GLDASpolyCentroid.Rd +++ b/man/GLDASpolyCentroid.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/GLDASpolyCentroid.R \name{GLDASpolyCentroid} \alias{GLDASpolyCentroid} -\title{Generate air temperature input files as well as air temperature stations file from NASA GLDAS remote sensing products.} +\title{Air temperature data from NASA GLDAS on centroid} \usage{ GLDASpolyCentroid( Dir = "./SWAT_INPUT/", @@ -15,9 +15,9 @@ GLDASpolyCentroid( \arguments{ \item{Dir}{A directory name to store gridded air temperature and air temperature stations files.} -\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'.} -\item{DEM}{A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{DEM}{A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'.} \item{start}{Beginning date for gridded air temperature data.} @@ -31,13 +31,13 @@ a scalar of maximum and minimum air temperature gridded data values at a pseudo This function downloads remote sensing data of \acronym{GLDAS} from \acronym{NASA} \acronym{GSFC} servers, extracts air temperature data from grids falling within a specified sub-basin(s) watershed shapefile, and assigns a pseudo air temperature gauge located at the centroid of the sub-basin(s) watershed a weighted-average daily minimum and maximum air temperature data. The function generates tables in a format that \acronym{SWAT} or other rainfall-runoff hydrological model requires for minimum and maximum air temperatures data input. The function also generates the air temperature stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that pseudo grids that correspond to the centroids of the watershed sub-basins. } \details{ -A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASpolyCentroid} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. +A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASpolyCentroid} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. The goal of the Global Land Data Assimilation System \acronym{GLDAS} is to ingest satellite and ground-based observational data products, using advanced land surface modeling and data assimilation techniques, in order to generate optimal fields of land surface states and fluxes (Rodell et al., 2004). \acronym{GLDAS} data set used in this function is the \acronym{GLDAS} Noah Land Surface Model L4 3 hourly 0.25 x 0.25 degree V2.1. The full suite of \acronym{GLDAS} data sets are available at \url{https://hydro1.gesdisc.eosdis.nasa.gov/dods/}. The \command{GLDASpolyCentroid} finds the minimum and maximum air temperatures for each day at each grid within the study watershed by searching for minima and maxima over the three hours air temperature data values available for each day and grid. -The \command{GLDASpolyCentroid} function relies on 'curl' tool to transfer data from \acronym{NASA} servers to a user machine, using HTTPS supported protocol. The 'curl' command embedded in this function to fetch \acronym{GLDAS} netcdf daily global files is designed to work seamlessly given that appropriate logging information are stored in the ".netrc" file and the cookies file ".urs_cookies" as explained in registering with the Earth Observing System Data and Information System. It is imperative to say here that a user machine should have 'curl' installed as a prerequisite to run \command{GLDASpolyCentroid}. +The \command{GLDASpolyCentroid} function relies on "curl" tool to transfer data from \acronym{NASA} servers to a user machine, using HTTPS supported protocol. The "curl" command embedded in this function to fetch \acronym{GLDAS} netcdf daily global files is designed to work seamlessly given that appropriate logging information are stored in the ".netrc" file and the cookies file ".urs_cookies" as explained in registering with the Earth Observing System Data and Information System. It is imperative to say here that a user machine should have "curl" installed as a prerequisite to run \command{GLDASpolyCentroid}. -The \acronym{GLDAS} V2.1 simulation started on January 1, 2000 using the conditions from the \acronym{GLDAS} V2.0 simulation. The \acronym{GLDAS} V2.1 simulation was forced with National Oceanic and Atmospheric Administration \acronym{NOAA}, Global Data Assimilation System \acronym{GDAS} atmospheric analysis fields (Derber et al., 1991), the disaggregated Global Precipitation Climatology Project \acronym{GPCP} precipitation fields (Adler et al., 2003), and the Air Force Weather Agency’s AGRicultural METeorological modeling system \acronym{AGRMET} radiation fields which became available for March 1, 2001 onward. +The \acronym{GLDAS} V2.1 simulation started on January 1, 2000 using the conditions from the \acronym{GLDAS} V2.0 simulation. The \acronym{GLDAS} V2.1 simulation was forced with National Oceanic and Atmospheric Administration \acronym{NOAA}, Global Data Assimilation System \acronym{GDAS} atmospheric analysis fields (Derber et al., 1991), the disaggregated Global Precipitation Climatology Project \acronym{GPCP} precipitation fields (Adler et al., 2003), and the Air Force Weather Agency's AGRicultural METeorological modeling system \acronym{AGRMET} radiation fields which became available for March 1, 2001 onward. } \note{ \command{start} should be equal to or greater than 2000-Jan-01. @@ -48,10 +48,10 @@ The \acronym{GLDAS} V2.1 simulation started on January 1, 2000 using the conditi DEM = "LowerMekong_dem.tif", start = "2015-12-1", end = "2015-12-3")} } \references{ -Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, \cr -doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. +Adler, R. F., G. J. Huffman, A. Chang, R. Ferraro, P.-P. Xie, J. Janowiak, B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. Arkin, and + E. Nelkin (2003), The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979–Present), J. Hydrometeorol., 4, 1147-1167, doi:10.1175/1525-7541(2003)004<1147:tvgpcp>2.0.co;2. -Adler, R. F., G. J. Huffman, A. Chang, R. Ferraro, P.-P. Xie, J. Janowiak, B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. Arkin, and E. Nelkin (2003), The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979–Present), J. Hydrometeorol., 4, 1147-1167, doi:10.1175/1525-7541(2003)004<1147:tvgpcp>2.0.co;2. +Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. Rodell, M., P. R. Houser, U. Jambor, J. Gottschalck, K. Mitchell, C.-J. Meng, K. Arsenault, B. Cosgrove, J. Radakovich, M. Bosilovich, J. K. Entin*, J. P. Walker, D. Lohmann, and D. Toll (2004), The Global Land Data Assimilation System, B. Am. Meteorol. Soc., 85, 381-394, doi:10.1175/bams-85-3-381. } diff --git a/man/GLDASwat.Rd b/man/GLDASwat.Rd index 4a4a8df..2a3c515 100644 --- a/man/GLDASwat.Rd +++ b/man/GLDASwat.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/GLDASwat.R \name{GLDASwat} \alias{GLDASwat} -\title{Generate SWAT air temperature input files as well as air temperature stations file from NASA GLDAS remote sensing products.} +\title{SWAT air temperature data from NASA GLDAS} \usage{ GLDASwat( Dir = "./SWAT_INPUT/", @@ -15,9 +15,9 @@ GLDASwat( \arguments{ \item{Dir}{A directory name to store gridded air temperature and air temperature stations files.} -\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'.} -\item{DEM}{A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{DEM}{A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'.} \item{start}{Beginning date for gridded air temperature data.} @@ -31,7 +31,7 @@ a scalar of maximum and minimum air temperature gridded data values at each poin This function downloads remote sensing data of \acronym{GLDAS} from \acronym{NASA} \acronym{GSFC} servers, extracts air temperature data from grids within a specified watershed shapefile, and then generates tables in a format that \acronym{SWAT} requires for minimum and maximum air temperature data input. The function also generates the air temperature stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that fall within the specified watershed. } \details{ -A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASwat} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. +A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{GLDAS} remote sensing data products at (\url{https://hydro1.gesdisc.eosdis.nasa.gov/data/GLDAS/GLDAS_NOAH025_3H.2.1/}). The function uses variable name ('Tair_f_inst') for air temperature in \acronym{GLDAS} data products. Units for gridded air temperature data are degrees in 'K'. The \command{GLDASwat} function outputs gridded air temperature (maximum and minimum) data in degrees 'C'. The goal of the Global Land Data Assimilation System \acronym{GLDAS} is to ingest satellite and ground-based observational data products, using advanced land surface modeling and data assimilation techniques, in order to generate optimal fields of land surface states and fluxes (Rodell et al., 2004). \acronym{GLDAS} dataset used in this function is the \acronym{GLDAS} Noah Land Surface Model L4 3 hourly 0.25 x 0.25 degree V2.1. The full suite of \acronym{GLDAS} datasets is available at \url{https://hydro1.gesdisc.eosdis.nasa.gov/dods/}. The \command{GLDASwat} finds the minimum and maximum air temperatures for each day at each grid within the study watershed by searching for minima and maxima over the three hours air temperature data values available for each day and grid. @@ -48,11 +48,10 @@ The \acronym{GLDAS} V2.1 simulation started on January 1, 2000 using the conditi DEM = "LowerMekong_dem.tif", start = "2015-12-1", end = "2015-12-3")} } \references{ -Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, \cr -doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. - Adler, R. F., G. J. Huffman, A. Chang, R. Ferraro, P.-P. Xie, J. Janowiak, B. Rudolf, U. Schneider, S. Curtis, D. Bolvin, A. Gruber, J. Susskind, P. Arkin, and E. Nelkin (2003), The Version-2 Global Precipitation Climatology Project (GPCP) Monthly Precipitation Analysis (1979–Present), J. Hydrometeorol., 4, 1147-1167, doi:10.1175/1525-7541(2003)004<1147:tvgpcp>2.0.co;2. +Derber, J. C., D. F. Parrish, and S. J. Lord (1991), The New Global Operational Analysis System at the National Meteorological Center, Weather Forecast, 6, 538-547, doi:10.1175/1520-0434(1991)006<0538:tngoas>2.0.co;2. + Rodell, M., P. R. Houser, U. Jambor, J. Gottschalck, K. Mitchell, C.-J. Meng, K. Arsenault, B. Cosgrove, J. Radakovich, M. Bosilovich, J. K. Entin*, J. P. Walker, D. Lohmann, and D. Toll (2004), The Global Land Data Assimilation System, B. Am. Meteorol. Soc., 85, 381-394, doi:10.1175/bams-85-3-381. } \author{ diff --git a/man/GPM_NRT.Rd b/man/GPM_NRT.Rd index fefcb51..03795e1 100644 --- a/man/GPM_NRT.Rd +++ b/man/GPM_NRT.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/GPM_NRT.R \name{GPM_NRT} \alias{GPM_NRT} -\title{Generate Near Real Time (NRT) rainfall from NASA GPM remote sensing products.} +\title{NASA GPM Near Real Time rainfall data} \usage{ GPM_NRT( Dir = "./INPUT/", @@ -15,9 +15,9 @@ GPM_NRT( \arguments{ \item{Dir}{A directory name to store gridded rainfall and rain stations files.} -\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'.} -\item{DEM}{A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{DEM}{A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'.} \item{start}{Beginning date for gridded rainfall data.} @@ -31,7 +31,7 @@ a scalar of rainfall gridded data values at each point within the study watershe This function downloads rainfall remote sensing data of \acronym{IMERG} from \acronym{NASA} \acronym{GSFC} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that any hydrological model requires for rainfall data input. The function also generates the rainfall stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that fall within the specified watershed. The minimum latency for this function is one day. } \details{ -A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDE.06/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products. Units for gridded rainfall data are 'mm'. +A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDE.06/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products. Units for gridded rainfall data are 'mm'. \acronym{IMERG} dataset is the GPM Level 3 \acronym{IMERG} *Early* Daily 0.1 x 0.1 deg (GPM_3IMERGDE) derived from the half-hourly \acronym{GPM_3IMERGHHE}. The derived result represents the final estimate of the daily accumulated precipitation. The dataset is produced at the \acronym{NASA} Goddard Earth Sciences (GES) Data and Information Services Center (DISC) by simply summing the valid precipitation retrievals for the day in GPM_3IMERGHHE and giving the result in (mm) \url{https://gpm.nasa.gov/data/directory}. diff --git a/man/GPMpolyCentroid.Rd b/man/GPMpolyCentroid.Rd index 11c9282..fc56de3 100644 --- a/man/GPMpolyCentroid.Rd +++ b/man/GPMpolyCentroid.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/GPMployCentroid.R \name{GPMpolyCentroid} \alias{GPMpolyCentroid} -\title{Generate rainfall input files as well as rain station file from NASA GPM remote sensing products.} +\title{NASA GPM rainfall data on centroid} \usage{ GPMpolyCentroid( Dir = "./SWAT_INPUT/", @@ -15,9 +15,9 @@ GPMpolyCentroid( \arguments{ \item{Dir}{A directory name to store gridded rainfall and rain stations files.} -\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'.} -\item{DEM}{A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{DEM}{A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'.} \item{start}{Beginning date for gridded rainfall data.} @@ -31,7 +31,7 @@ a scalar of rainfall gridded data values at a pseudo rain grid located at the ce This function downloads rainfall remote sensing data of \acronym{TRMM} and \acronym{IMERG} from \acronym{NASA} \acronym{GSFC} servers, extracts data from grids falling within a specified sub-basin(s) watershed shapefile and assigns a pseudo rainfall gauge located at the centroid of the sub-basin(s) watershed a weighted-average daily rainfall data. The function generates rainfall tables in a format that \acronym{SWAT} or other rainfall-runoff hydrological model requires for rainfall data input. The function also generates the rainfall stations file summary input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those pseudo grids that correspond to the centroids of the watershed sub-basins. } \details{ -A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. +A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. \acronym{IMERG} dataset is the GPM Level 3 \acronym{IMERG} *Final* Daily 0.1 x 0.1 deg (GPM_3IMERGDF) derived from the half-hourly GPM_3IMERGHH. The derived result represents the final estimate of the daily accumulated precipitation. The dataset is produced at the \acronym{NASA} Goddard Earth Sciences (GES) Data and Information Services Center (DISC) by simply summing the valid precipitation retrievals for the day in GPM_3IMERGHH and giving the result in (mm) \url{https://gpm.nasa.gov/data/directory}. diff --git a/man/GPMswat.Rd b/man/GPMswat.Rd index 71d0799..633e41c 100644 --- a/man/GPMswat.Rd +++ b/man/GPMswat.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/GPMswat.R \name{GPMswat} \alias{GPMswat} -\title{Generate SWAT rainfall input files as well as rain stations file from NASA GPM remote sensing products.} +\title{SWAT rainfall data from NASA GPM} \usage{ GPMswat( Dir = "./SWAT_INPUT/", @@ -15,9 +15,9 @@ GPMswat( \arguments{ \item{Dir}{A directory name to store gridded rainfall and rain stations files.} -\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection crs='+proj=longlat +datum=WGS84'.} -\item{DEM}{A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +datum=WGS84').} +\item{DEM}{A study watershed digital elevation model raster in a geographic projection crs='+proj=longlat +datum=WGS84'.} \item{start}{Beginning date for gridded rainfall data.} @@ -31,7 +31,7 @@ a scalar of rainfall gridded data values at each point within the study watershe This function downloads rainfall remote sensing data of \acronym{TRMM} and \acronym{IMERG} from \acronym{NASA} \acronym{GSFC} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that \acronym{SWAT} requires for rainfall data input. The function also generates the rainfall stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected grids that fall within the specified watershed. } \details{ -A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. +A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} GESDISC Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server address for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server address for \acronym{TRMM} remote sensing data products (\url{https://disc2.gesdisc.eosdis.nasa.gov/data/TRMM_RT/TRMM_3B42RT_Daily.7/}). The function uses variable name ('precipitationCal') for rainfall in \acronym{IMERG} data products and variable name ('precipitation') for \acronym{TRMM} rainfall data products. Units for gridded rainfall data are 'mm'. \acronym{IMERG} dataset is the GPM Level 3 \acronym{IMERG} *Final* Daily 0.1 x 0.1 deg (GPM_3IMERGDF) derived from the half-hourly GPM_3IMERGHH. The derived result represents the final estimate of the daily accumulated precipitation. The dataset is produced at the \acronym{NASA} Goddard Earth Sciences (GES) Data and Information Services Center (DISC) by simply summing the valid precipitation retrievals for the day in GPM_3IMERGHH and giving the result in (mm) \url{https://gpm.nasa.gov/data/directory}. diff --git a/man/NEX_GDDP_CMIP5.Rd b/man/NEX_GDDP_CMIP5.Rd index c7aa474..36feccc 100644 --- a/man/NEX_GDDP_CMIP5.Rd +++ b/man/NEX_GDDP_CMIP5.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/NEXGDDP.R \name{NEX_GDDP_CMIP5} \alias{NEX_GDDP_CMIP5} -\title{Generate rainfall or air temperature as well as climate input stations file from NASA NEX-GDDP remote sensing climate change data products needed to drive various hydrological models.} +\title{CMIP5 climate data from NASA NEX-GDDP} \usage{ NEX_GDDP_CMIP5( Dir = "./INPUT/", @@ -18,9 +18,9 @@ NEX_GDDP_CMIP5( \arguments{ \item{Dir}{A directory name to store gridded climate data and stations files.} -\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs').} +\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection crs ='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'.} -\item{DEM}{A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs').} +\item{DEM}{A study watershed digital elevation model raster in a geographic projection crs ='+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'.} \item{start}{Beginning date for gridded climate data.} @@ -40,7 +40,7 @@ a scalar of climate change gridded data values at each point within the study wa This function downloads climate change data of rainfall and air temperature from \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP} \acronym{GSFC} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that any hydrological model requires for rainfall or air temperature data input. The function also generates the climate stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected climatological grids that fall within the specified watershed. The \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP} dataset is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 5 \acronym{CMIP5} and across two of the four greenhouse gas emissions scenarios known as Representative Concentration Pathways \acronym{RCPs} (rcp45, rcp85). } \details{ -A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server for \acronym{NEX-GDDP} climate change data products at (\url{https://www.nccs.nasa.gov/services/climate-data-services}). The function uses variable name ('pr') for rainfall in \acronym{NEX-GDDP} data products and variable name ('tas') for \acronym{NEX-GDDP} minimum ('tasmin') and maximum ('tasmax') air temperature data products. The \command{NEX-GDDP} function outputs gridded rainfall data in 'mm' and gridded air temperature (maximum and minimum) data in degrees 'C'. +A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} Goddard Space Flight Center server for \acronym{NEX-GDDP} climate change data products at (\url{https://www.nccs.nasa.gov/services/climate-data-services}). The function uses variable name ('pr') for rainfall in \acronym{NEX-GDDP} data products and variable name ('tas') for \acronym{NEX-GDDP} minimum ('tasmin') and maximum ('tasmax') air temperature data products. The \command{NEX-GDDP} function outputs gridded rainfall data in 'mm' and gridded air temperature (maximum and minimum) data in degrees 'C'. \acronym{NEX-GDDP} dataset is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 5 \acronym{CMIP5} (Taylor et al. 2012) and across two of the four greenhouse gas emissions scenarios known as Representative Concentration Pathways \acronym{RCPs} (Meinshausen et al. 2011). The \acronym{CMIP5} \acronym{GCM} runs were developed in support of the Fifth Assessment Report of the Intergovernmental Panel on Climate Change \acronym{IPCC AR5}. This dataset includes downscaled projections from the 21 models and scenarios for which daily scenarios were produced and distributed under \acronym{CMIP5}. The Bias-Correction Spatial Disaggregation \acronym{BCSD} method used in generating the \acronym{NEX-GDDP} dataset is a statistical downscaling algorithm specifically developed to address the current limitations of the global \acronym{GCM} outputs (Wood et al. 2002; Wood et al. 2004; Maurer et al. 2008; Thrasher et al. 2012). The \acronym{NEX-GDDP} climate projections is downscaled at a spatial resolution of 0.25 degrees x 0.25 degrees (approximately 25 km x 25 km). The \command{NEX_GDDP_CMIP5} downscales the \acronym{NEX-GDDP} data to grid points of 0.1 degrees x 0.1 degrees following nearest point methods described by Mohammed et al. (2018). diff --git a/man/NEX_GDDP_CMIP6.Rd b/man/NEX_GDDP_CMIP6.Rd index 11484fb..f7c3908 100644 --- a/man/NEX_GDDP_CMIP6.Rd +++ b/man/NEX_GDDP_CMIP6.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/NEXGDDP_CMIP6.R \name{NEX_GDDP_CMIP6} \alias{NEX_GDDP_CMIP6} -\title{Generate rainfall or air temperature as well as climate input stations file from NASA NEX-GDDP-CMIP6 remote sensing climate change data products needed to drive various hydrological models.} +\title{CMIP6 climate data from NASA NEX-GDDP} \usage{ NEX_GDDP_CMIP6( Dir = "./INPUT/", @@ -18,9 +18,9 @@ NEX_GDDP_CMIP6( \arguments{ \item{Dir}{A directory name to store gridded climate data and stations files.} -\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs').} +\item{watershed}{A study watershed shapefile spatially describing polygon(s) in a geographic projection crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'.} -\item{DEM}{A study watershed digital elevation model raster in a geographic projection sp::CRS('+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs').} +\item{DEM}{A study watershed digital elevation model raster in a geographic projection crs = '+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs'.} \item{start}{Beginning date for gridded climate data.} @@ -40,14 +40,14 @@ a scalar of climate change gridded data values at each point within the study wa This function downloads climate change data of rainfall and air temperature from \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP-CMIP6} \acronym{AMES} servers, extracts data from grids within a specified watershed shapefile, and then generates tables in a format that any hydrological model requires for rainfall or air temperature data input. The function also generates the climate stations file input (file with columns: ID, File NAME, LAT, LONG, and ELEVATION) for those selected climatological grids that fall within the specified watershed. The \acronym{NASA} Earth Exchange Global Daily Downscaled Projections \acronym{NEX-GDDP-CMIP6} data set is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 6 \acronym{CMIP6} and across two of the four "Tier 1" greenhouse gas emissions scenarios known as Shared Socioeconomic Pathways \acronym{SSPs}. } \details{ -A user should visit \url{https://disc.gsfc.nasa.gov/data-access} to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. +A user should visit \url{https://disc.gsfc.nasa.gov/information/documents} Data Access document to register with the Earth Observing System Data and Information System (\acronym{NASA Earthdata}) and then authorize \acronym{NASA} \acronym{GESDISC} Data Access to successfully work with this function. The function accesses \acronym{NASA} Goddard Space Flight Center server for \acronym{IMERG} remote sensing data products at (\url{https://gpm1.gesdisc.eosdis.nasa.gov/data/GPM_L3/GPM_3IMERGDF.06/}), and \acronym{NASA} AMES Research Center server for \acronym{NEX-GDDP-CMIP6} climate change data products at \url{https://www.nccs.nasa.gov/services/data-collections/land-based-products/nex-gddp-cmip6}. The function uses variable name ('pr') for rainfall in \acronym{NEX-GDDP-CMIP6} data products and variable name ('tas') for \acronym{NEX-GDDP-CMIP6} minimum ('tasmin') and maximum ('tasmax') air temperature data products. The \command{NEX_GDDP_CMIP6} function outputs gridded rainfall data in 'mm' and gridded air temperature (maximum and minimum) data in degrees 'C'. \acronym{NEX-GDDP-CMIP6} dataset is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model \acronym{GCM} runs conducted under the Coupled Model Intercomparison Project Phase 6 \acronym{CMIP6} (Eyring et al. 2016) and across the four "Tier 1" greenhouse gas emissions scenarios known as Shared Socioeconomic Pathways \acronym{SSPs} (O'Neil et al. 2016; Meinshausen et al. 2020). The \acronym{CMIP6} \acronym{GCM} runs were developed in support of the Sixth Assessment Report -of the Intergovernmental Panel on Climate Change \acronym{IPCC AR6}. This data set includes downscaled projections from the 35 models and scenarios for which daily scenarios were produced and distributed under \acronym{CMIP6}. Please visit \acronym{NCCS} Dataportal - Datashare \url{https://ds.nccs.nasa.gov/thredds2/catalog/AMES/NEX/GDDP-CMIP6/catalog.html} to ensure that the requested model and greenhouse gas emissions scenario (\acronym{SSPs}) is available on the server before you execute your run. +of the Intergovernmental Panel on Climate Change \acronym{IPCC AR6}. This data set includes downscaled projections from the 35 models and scenarios for which daily scenarios were produced and distributed under \acronym{CMIP6}. Please visit \acronym{NCCS} Dataportal - Datashare technical note \url{https://www.nccs.nasa.gov/sites/default/files/NEX-GDDP-CMIP6-Tech_Note.pdf} to ensure that the requested model and greenhouse gas emissions scenario (\acronym{SSPs}) is available on the server before you execute your run. The Bias-Correction Spatial Disaggregation \acronym{BCSD} method used in generating the \acronym{NEX-GDDP-CMIP6} data set is a statistical downscaling algorithm specifically developed to address the current limitations of the global \acronym{GCM} outputs (Wood et al. 2002; Wood et al. 2004; Maurer et al. 2008; Thrasher et al. 2012). The \acronym{NEX-GDDP-CMIP6} climate projections is downscaled at a spatial resolution of 0.25 degrees x 0.25 degrees (approximately 25 km x 25 km). The \command{NEX_GDDP_CMIP6} downscales the \acronym{NEX-GDDP-CMIP6} data to grid points of 0.1 degrees x 0.1 degrees following nearest point methods described by Mohammed et al. (2018). diff --git a/vignettes/About.Rmd b/vignettes/About.Rmd index aed73ff..d75938c 100644 --- a/vignettes/About.Rmd +++ b/vignettes/About.Rmd @@ -17,10 +17,13 @@ knitr::opts_chunk$set( ``` The _NASAaccess_ package is a simple and needed tool for accessing and processing remote sensing data. -The _NASAaccess_ package has multiple functions that help the user to access and reformat hydrological data for easy ingest into various hydrological models. Since the package functions touch [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") data repositories to retrieve data, the user must set up a registration account with [Earthdata](https://urs.earthdata.nasa.gov/users/new) as well as authorizing [NASA GES DISC](https://disc.gsfc.nasa.gov/data-access) data access. The package user should make sure that his(her) local machines has [curl](https://curl.se/) installed properly. Further instructions on creating the ".netrc" and ".urs_cookies" files can be accessed at [*How To Access Data NASA data With cURL And Wget wiki page*](https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+cURL+And+Wget). Creating the ".netrc" file at the user machine 'Home' directory and storing the [NASA GES DISC](https://disc.gsfc.nasa.gov/) user logging information in it is needed to execute the package commands. +The _NASAaccess_ package has multiple functions that help the user to access and reformat hydrological data for easy ingest into various hydrological models. Since the package functions touch [NASA](https://www.nasa.gov/ "The National Aeronautics and Space Administration") data repositories to retrieve data, the user must set up a registration account with [Earthdata](https://urs.earthdata.nasa.gov/users/new) as well as authorizing [NASA GES DISC](https://disc.gsfc.nasa.gov/information/documents?title=Data%20Access) data access. The package user should make sure that his(her) local machines has [curl](https://curl.se/) installed properly. Further instructions on creating the ".netrc" and ".urs_cookies" files can be accessed at [*How To Access Data NASA data With cURL And Wget wiki page*](https://wiki.earthdata.nasa.gov/display/EL/How+To+Access+Data+With+cURL+And+Wget). Creating the ".netrc" file at the user machine 'Home' directory and storing the [NASA GES DISC](https://disc.gsfc.nasa.gov/) user logging information in it is needed to execute the package commands. Note: for Windows users the [NASA GES DISC](https://disc.gsfc.nasa.gov/) logging information should be saved in a file named "_netrc" beside the ".netrc" one. +## Video Abstract + +[Video abstract for "Technical note: NASAaccess – A tool for access, reformatting, and visualization of remotely sensed earth observation and climate data"](https://doi.org/10.5446/63008). ## Index @@ -30,9 +33,15 @@ Note: for Windows users the [NASA GES DISC](https://disc.gsfc.nasa.gov/) logging * [CMIP6 Climate](NEXGDDP-CMIP6.html) + ## Built with ```{r} sessionInfo() ``` + + +## Reference + +Mohammed, I.N., Bustamante, E.G.R., Bolten, J.D., Nelson, E.J., 2023. Technical note: NASAaccess – a tool for access, reformatting, and visualization of remotely sensed earth observation and climate data. Hydrol. Earth Syst. Sci. 27, 3621-3642, https://doi.org/10.5194/hess-27-3621-2023 diff --git a/vignettes/GLDAS.Rmd b/vignettes/GLDAS.Rmd index c39963e..19a4ee8 100644 --- a/vignettes/GLDAS.Rmd +++ b/vignettes/GLDAS.Rmd @@ -19,18 +19,19 @@ knitr::opts_chunk$set( ) - -library(raster) -library(sf) - #Reading input data -dem_path <- system.file("extdata", "DEM_TX.tif", package = "NASAaccess") -shape_path <- system.file("extdata", "basin.shp", package = "NASAaccess") +dem_path <- system.file("extdata", + "DEM_TX.tif", + package = "NASAaccess") +shape_path <- system.file("extdata", + "basin.shp", + package = "NASAaccess") +library(terra) -dem <- raster::raster(dem_path) -shape <- sf::st_read(shape_path) +dem <- terra::rast(dem_path) +shape <- terra::vect(shape_path) ``` @@ -43,9 +44,19 @@ Let's explore `GLDASpolyCentroid` and `GLDASswat` functions. ## Basic use -Let's use the example watersheds that we introduced with `GPMswat` and `GPMpolyCentroid`. Please visit _NASAaccess_ [GPM](.\GPM.pdf) functions for more information. +Let's use the example watersheds that we introduced with `GPMswat` and `GPMpolyCentroid`. Please visit _NASAaccess_ [GPM](GPM.html) functions for more information. ```{r, eval=FALSE, echo=TRUE} +#Reading input data + +dem_path <- system.file("extdata", + "DEM_TX.tif", + package = "NASAaccess") + +shape_path <- system.file("extdata", + "basin.shp", + package = "NASAaccess") + library(NASAaccess) GLDASwat(Dir = "./GLDASwat/", @@ -76,10 +87,11 @@ Now, let's see the location of the `GLDASwat` generated grid points ```{r} library(ggplot2) -ggplot() + - geom_sf(data = shape, - fill = NA, - colour = 'red') + +library(tidyterra) +ggplot(shape) + + + geom_spatvector(color='red',fill=NA) + + geom_point(data=GLDASwat.table,aes(x=LONG,y=LAT)) @@ -124,6 +136,7 @@ Now let's examine the `GLDASpolyCentroid` generated outputs ```{r} library(ggplot2) +library(tidyterra) GLDASpolyCentroid.tempMaster <- system.file('extdata/GLDASpolyCentroid', 'temp_Master.txt', package = 'NASAaccess') @@ -131,10 +144,10 @@ GLDASpolyCentroid.tempMaster <- system.file('extdata/GLDASpolyCentroid', GLDASpolyCentroid.table<-read.csv(GLDASpolyCentroid.tempMaster) #plot -ggplot() + - geom_sf(data = shape, - fill = NA, - colour = 'red') + +ggplot(shape) + + + geom_spatvector(color='red',fill=NA) + + geom_point(data=GLDASpolyCentroid.table, aes(x=LONG,y=LAT)) + coord_sf(xlim=c(-96,-95.2),ylim=c(29.7,30)) diff --git a/vignettes/GPM.Rmd b/vignettes/GPM.Rmd index 22503ce..0669554 100644 --- a/vignettes/GPM.Rmd +++ b/vignettes/GPM.Rmd @@ -27,11 +27,9 @@ Let's explore `GPMpolyCentroid` and `GPMswat` functions. Let's look at an example watershed that we want to examine near Houston, TX: ```{r} -library(ggmap) -library(raster) +library(leaflet) library(ggplot2) -library(rgdal) - +library(terra) #Reading input data dem_path <- system.file("extdata", "DEM_TX.tif", @@ -42,50 +40,35 @@ shape_path <- system.file("extdata", package = "NASAaccess") -dem <- raster::raster(dem_path) +dem <- terra::rast(dem_path) + +shape <- terra::vect(shape_path) -shape <- rgdal::readOGR(shape_path) -shape.df <- ggplot2::fortify(shape) #plot the watershed data -myMap <- get_stamenmap(bbox = c(left = -96, - bottom = 29.7, - right = -95.2, - top = 30), - maptype = "terrain", - crop = TRUE, - zoom = 10) +myMap <- leaflet() %>% + addTiles() %>% + fitBounds(-96, 29.7, -95.2, 30) %>% + addPolygons(lng=terra::crds(shape)[,1], + lat=terra::crds(shape)[,2]) -ggmap(myMap) + - geom_polygon(data = shape.df, - aes(x = long, y = lat, group = group), - fill = NA, size = 0.5, color = 'red') +myMap ``` The geographic layout of the White Oak Bayou watershed example used in this demonstration is depicted above. Whiteoak Bayou is a tributary for the Buffalo Bayou River (Harris County, Texas). In order to use *NASAaccess* we also need a digital elevation model (DEM) raster layer. Let's see the White Oak Bayou watershed DEM and a more closer look at the study watershed example. ```{r} -# create a plot of our DEM raster along with watershed - -plot(dem, - main="White Oak Bayou Watershed with Digital Elevation Model (DEM)", - col=rev(bpy.colors()), - xlab='lon', - ylab='lat', - legend = T, - legend.args=list(text='Elevation (m)', - side=4, - font=2, - line=2.5, - cex=0.8)) +# create a plot of our DEM raster along with watershed (i.e., elevation in meters) + +plot(dem,main="White Oak Bayou Watershed with Digital Elevation Model (DEM)") plot(shape , add = TRUE) @@ -125,11 +108,9 @@ Now, let's see the location of these generated grid points: ```{r} library(ggplot2) -ggplot() + - geom_polygon(data = shape.df, - aes(x = long, y = lat, group = group), - fill = NA, - colour = 'red') + +library(tidyterra) +ggplot(shape) + + geom_spatvector(color='red',fill=NA) + geom_point(data=GPMswat.table, aes(x=LONG, y=LAT, @@ -169,6 +150,7 @@ Examining the rainfall station file generated by `GPMpolyCentroid` ```{r} library(ggplot2) +library(tidyterra) GPMpolyCentroid.precipitationMaster <- system.file('extdata/GPMpolyCentroid', 'precipitationMaster.txt', package = 'NASAaccess') @@ -177,11 +159,10 @@ GPMpolyCentroid.precipitation.table <- read.csv(GPMpolyCentroid.precipitationMas # plotting -ggplot() + - geom_polygon(data = shape.df, - aes(x = long, y = lat, group = group), - fill = NA, - colour = 'red') + +ggplot(shape) + + + geom_spatvector(color='red',fill=NA) + + geom_point(data=GPMpolyCentroid.precipitation.table, aes(x=LONG,y=LAT)) ``` @@ -201,7 +182,8 @@ GPMpolyCentroid.precipitation.data <- read.csv(GPMpolyCentroid.precipitation.rec # since data started on 2019-08-01 days <- seq.Date(from = as.Date('2019-08-01'), - length.out = dim(GPMpolyCentroid.precipitation.data)[1], by = 'day') + length.out = dim(GPMpolyCentroid.precipitation.data)[1], + by = 'day') # plotting the precipitation time series diff --git a/vignettes/NEXGDDP-CMIP6.Rmd b/vignettes/NEXGDDP-CMIP6.Rmd index 6b7ba0e..57139fc 100644 --- a/vignettes/NEXGDDP-CMIP6.Rmd +++ b/vignettes/NEXGDDP-CMIP6.Rmd @@ -16,22 +16,8 @@ knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) - -library(raster) -library(sf) - -#Reading input data -dem_path <- system.file("extdata", "DEM_TX.tif", package = "NASAaccess") -shape_path <- system.file("extdata", "basin.shp", package = "NASAaccess") - - -dem <- raster::raster(dem_path) -shape <- sf::st_read(shape_path) ``` - - - _NASAaccess_ has a handy tool to access, extract, and reformat climate change data of rainfall and air temperature from [NASA Earth Exchange Global Daily Downscaled Projections NEX-GDDP-CMIP6 AMES servers](https://www.nccs.nasa.gov/services/data-collections/land-based-products/nex-gddp-cmip6) for grids within a specified watershed. [NEX-GDDP-CMIP6](https://www.nccs.nasa.gov/sites/default/files/NEX-GDDP-CMIP6-Tech_Note.pdf) dataset is comprised of downscaled climate scenarios for the globe that are derived from the General Circulation Model GCM runs conducted under the Coupled Model Intercomparison Project Phase 6 CMIP6 [@RN1450] and across two of the four "Tier 1" greenhouse gas emissions scenarios known as Shared Socioeconomic Pathways (SSPs) [@RN1451; @RN1452]. The CMIP6 GCM runs were developed in support of the Sixth Assessment Report of the Intergovernmental Panel on Climate Change IPCC AR6. This dataset includes downscaled projections from the [35 models and scenarios](https://pcmdi.llnl.gov/CMIP6/) for which daily scenarios were produced and distributed under CMIP6. @@ -40,19 +26,33 @@ The Bias-Correction Spatial Disaggregation BCSD method used in generating the [N ## Basic use -Let's use the example watersheds that we introduced with `GPMswat` and `GPMpolyCentroid`. Please visit _NASAaccess_ [GPM](.\GPM.pdf) functions for more information. +Let's use the example watersheds that we introduced with `GPMswat` and `GPMpolyCentroid`. Please visit _NASAaccess_ [GPM](GPM.html) functions for more information. + ```{r, eval=FALSE, echo=TRUE} +#Reading input data + +dem_path <- system.file("extdata", + "DEM_TX.tif", + package = "NASAaccess") + +shape_path <- system.file("extdata", + "basin.shp", + package = "NASAaccess") + +#CMIP6 example for air temperature + library(NASAaccess) -NEX_GDDP_CMIP6(Dir = "./NEX_GDDP_CMIP6/", +NEX_GDDP_CMIP6( + Dir = "./NEX_GDDP_CMIP6/", watershed = shape_path, - dem_path, - start = "2060-12-1", - end = "2060-12-3", - model = 'ACCESS-CM2', - type = 'tas', - slice = 'ssp245') + DEM = dem_path, + start = "2060-12-1", + end = "2060-12-3", + model = 'ACCESS-CM2', + type = 'tas', + slice = 'ssp245') ``` Let's examine the air temperature station file diff --git a/vignettes/NEXGDDP.Rmd b/vignettes/NEXGDDP.Rmd index 058a737..15c3acb 100644 --- a/vignettes/NEXGDDP.Rmd +++ b/vignettes/NEXGDDP.Rmd @@ -18,17 +18,6 @@ knitr::opts_chunk$set( fig.height = 7, comment = "#>" ) - -library(raster) -library(sf) - -#Reading input data -dem_path <- system.file("extdata", "DEM_TX.tif", package = "NASAaccess") -shape_path <- system.file("extdata", "basin.shp", package = "NASAaccess") - - -dem <- raster::raster(dem_path) -shape <- sf::st_read(shape_path) ``` @@ -40,19 +29,29 @@ The Bias-Correction Spatial Disaggregation BCSD method used in generating the [N ## Basic use -Let's use the example watersheds that we introduced with `GPMswat` and `GPMpolyCentroid`. Please visit _NASAaccess_ [GPM](.\GPM.pdf) functions for more information. +Let's use the example watersheds that we introduced with `GPMswat` and `GPMpolyCentroid`. Please visit _NASAaccess_ [GPM](GPM.html) functions for more information. ```{r, eval=FALSE, echo=TRUE} +#Reading input data +dem_path <- system.file("extdata", + "DEM_TX.tif", + package = "NASAaccess") + +shape_path <- system.file("extdata", + "basin.shp", + package = "NASAaccess") + + library(NASAaccess) NEX_GDDP_CMIP5(Dir = "./NEX_GDDP_CMIP5/", - watershed = shape_path, - dem_path, + watershed = shape_path, + DEM = dem_path, start = "2060-12-1", end = "2060-12-3", model = 'IPSL-CM5A-MR', - type = 'pr', - slice = 'rcp85') + type = 'pr', + slice = 'rcp85') ``` Let's examine the precipitation station file