In this example we want to download, import and plot over the 3500 km2 wide region of interest (mapped below) :

  • a 30 days-long time series of land surface temperature from MODIS Terra LST (spatial resolution : 1 km ; temporal resolution : 1 day),
  • the same 30 days-long times series of precipitations from Global Precipitation Measurement (GPM) (spatial resolution : 1° ; temporal resolution : 1 day)
  • the same 30 days-long times series of soil moisture from SMAP Daily (spatial resolution : 9 km ; temporal resolution : 2/3 days)

Check which collections are available

Before starting, let’s identify our collections of interest among the collections that are available for download, using the function odr_list_collections() :

collections <- opendapr::odr_list_collections()
print(str(collections))
#> 'data.frame':    77 obs. of  26 variables:
#>  $ collection              : chr  "GPM_3IMERGDE.06" "GPM_3IMERGDF.06" "GPM_3IMERGDL.06" "GPM_3IMERGHH.06" ...
#>  $ source                  : chr  "GPM" "GPM" "GPM" "GPM" ...
#>  $ long_name               : chr  "GPM IMERG Early Precipitation L3 1 day 0.1 degree x 0.1 degree V06" "GPM IMERG Final Precipitation L3 1 day 0.1 degree x 0.1 degree V06" "GPM IMERG Late Precipitation L3 1 day 0.1 degree x 0.1 degree V06" "GPM IMERG Final Precipitation L3 Half Hourly 0.1 degree x 0.1 degree V06" ...
#>  $ nature                  : chr  "Rainfall" "Rainfall" "Rainfall" "Rainfall" ...
#>  $ provider                : chr  "NASA GES DISC" "NASA GES DISC" "NASA GES DISC" "NASA GES DISC" ...
#>  $ url_metadata            : chr  "https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGDE_06/summary" "https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGDF_06/summary" "https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGDL_06/summary" "https://disc.gsfc.nasa.gov/datasets/GPM_3IMERGHH_06/summary" ...
#>  $ doi                     : chr  "https://doi.org/10.5067/GPM/IMERGDE/DAY/06" "https://doi.org/10.5067/GPM/IMERGDF/DAY/06" "https://doi.org/10.5067/GPM/IMERGDL/DAY/06" "https://doi.org/10.5067/GPM/IMERG/3B-HH/06" ...
#>  $ version                 : num  6 6 6 6 6 6 6 6 6 6 ...
#>  $ spatial_resolution_m    : int  10000 10000 10000 10000 10000 10000 10000 500 500 500 ...
#>  $ temporal_resolution     : int  1 1 1 30 30 30 1 365 8 4 ...
#>  $ temporal_resolution_unit: chr  "day" "day" "day" "minute" ...
#>  $ spatial_coverage        : chr  "Global" "Global" "Global" "Global" ...
#>  $ start_date              : chr  "2000-06-01" "2000-06-01" "2000-06-01" "2000-06-01" ...
#>  $ end_date                : chr  "ongoing" "ongoing" "ongoing" "ongoing" ...
#>  $ indicative_latency_days : int  1 100 1 100 0 1 100 774 12 8 ...
#>  $ url_manual_access       : chr  "https://search.earthdata.nasa.gov/search?q=GPM_3IMERGDE_06" "https://search.earthdata.nasa.gov/search?q=GPM_3IMERGDF_06" "https://search.earthdata.nasa.gov/search?q=GPM_3IMERGDL_06" "https://search.earthdata.nasa.gov/search?q=GPM_3IMERGHH_06" ...
#>  $ status                  : chr  "Implemented" "Implemented" "Implemented" "Implemented" ...
#>  $ login                   : chr  "earthdata" "earthdata" "earthdata" "earthdata" ...
#>  $ url_opendapserver       : chr  "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/" "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/" "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/" "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/" ...
#>  $ url_programmatic_access : chr  "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDE.06/" "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDF.06/" "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDL.06/" "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGHH.06/" ...
#>  $ url_opendapexample      : chr  "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDE.06/2000/06/3B-DAY-E.MS.MRG.3IMERG.20000601-S00"| __truncated__ "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/hyrax/GPM_L3/GPM_3IMERGDF.06/2000/06/3B-DAY.MS.MRG.3IMERG.20000601"| __truncated__ "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/hyrax/GPM_L3/GPM_3IMERGDL.06/2000/06/3B-DAY-L.MS.MRG.3IMERG.200006"| __truncated__ "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/hyrax/GPM_L3/GPM_3IMERGHH.06/2000/153/3B-HHR.MS.MRG.3IMERG.2000060"| __truncated__ ...
#>  $ dim_lon                 : chr  "lon" "lon" "lon" "lon" ...
#>  $ dim_lat                 : chr  "lat" "lat" "lat" "lat" ...
#>  $ dim_time                : chr  "time" "time" "time" "time" ...
#>  $ dim_proj                : chr  NA NA NA NA ...
#>  $ crs                     : chr  "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,3" "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0" "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,4" ...
#> NULL

So our collections of interest are :

  • MOD11A1.006
  • GPM_3IMERGDF.06
  • SPL3SMP_E.003

Setup the region and time range of interest

First we call the useful packages and we define the ROI and the time frame.

require(opendapr)
require(sf)
require(stars)
require(raster)
require(ncdf4)
require(magrittr)
require(purrr)
# Set ROI and time range of interest
roi <- st_as_sf(data.frame(geom="POLYGON ((-5.82 9.54, -5.42 9.55, -5.41 8.84, -5.81 8.84, -5.82 9.54))"), wkt="geom", crs = 4326)
time_range <- as.Date(c("2017-01-01","2017-01-30"))

Login to EOSDIS Earthdata

And we login to EOSDIS Earthdata with our credentials. To create an account go to : https://urs.earthdata.nasa.gov/.

Get the OPeNDAP URLs of the data to download

With the function odr_get_url(), we get the https URLs for our collections of interest (MOD11A1.006, GPM_3IMERGDF.06 and SPL3SMP_E.003) given our ROI and time frame.

Note on the use of the parameter variables of the function odr_get_url() :

The collections generally contain several variables (sometimes called “dimensions”, or “bands”). As an example, MODIS LST products contain one band for the day temperature (“LST_Day_1km”), one band for the night temperature (“LST_Night_1km”), etc. To get all the variables available for a given collection along with information for each of them (description, etc.), use the function odr_list_variables(). e.g. :

tail(odr_list_variables("MOD11A1.006"))
#>             name                                             long_name
#> 13 Clear_day_cov                                day clear-sky coverage
#> 14 Day_view_angl View zenith angle of daytime Land-surface Temperature
#> 15   LST_Day_1km       Daily daytime 1km grid Land-surface Temperature
#> 16       Emis_31                                    Band 31 emissivity
#> 17       Emis_32                                    Band 32 emissivity
#> 18          time                                                  <NA>
#>    units                                   indices
#> 13  <NA> [time = 6845] [YDim = 1200] [XDim = 1200]
#> 14   deg [time = 6845] [YDim = 1200] [XDim = 1200]
#> 15     K [time = 6845] [YDim = 1200] [XDim = 1200]
#> 16  <NA> [time = 6845] [YDim = 1200] [XDim = 1200]
#> 17  <NA> [time = 6845] [YDim = 1200] [XDim = 1200]
#> 18  <NA>                             [time = 6845]
#>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                      all_info
#> 13                                                                           Array of 32 bit Reals [time = 0..6844][YDim = 0..1199][XDim = 0..1199]coordinates: Latitude Longitude\nlong_name: day clear-sky coverage\nday_clear_sky_cov: day_clear_sky_cov data * scale_factor\nscale_factor_err: 0.0000000000000000\nadd_offset_err: 0.0000000000000000\ncalibrated_nt: 5\norig_scale_factor: 0.00050000000000000001\n_FillValue: 0\nvalid_min: 0.0005000000237\nvalid_max: 32.76750183\nNumber_Type_Orig: uint16\ngrid_mapping: MODIS_Grid_Daily_1km_LST_eos_cf_projection
#> 14 Array of 32 bit Reals [time = 0..6844][YDim = 0..1199][XDim = 0..1199]coordinates: Latitude Longitude\nlong_name: View zenith angle of daytime Land-surface Temperature\nunits: deg\nView_angle: View_angle data * scale_factor + add_offset\nscale_factor_err: 0.0000000000000000\nadd_offset_err: 0.0000000000000000\ncalibrated_nt: 5\norig_scale_factor: 1.0000000000000000\norig_add_offset: -65.000000000000000\n_FillValue: 255\nvalid_min: -65.00000000\nvalid_max: 65.00000000\nNumber_Type_Orig: uint8\ngrid_mapping: MODIS_Grid_Daily_1km_LST_eos_cf_projection
#> 15                                                                          Array of 32 bit Reals [time = 0..6844][YDim = 0..1199][XDim = 0..1199]coordinates: Latitude Longitude\nlong_name: Daily daytime 1km grid Land-surface Temperature\nunits: K\nLST: LST data * scale_factor\nscale_factor_err: 0.0000000000000000\nadd_offset_err: 0.0000000000000000\ncalibrated_nt: 5\norig_scale_factor: 0.020000000000000000\n_FillValue: 0\nvalid_min: 150.0000000\nvalid_max: 1310.699951\nNumber_Type_Orig: uint16\ngrid_mapping: MODIS_Grid_Daily_1km_LST_eos_cf_projection
#> 16                                                     Array of 32 bit Reals [time = 0..6844][YDim = 0..1199][XDim = 0..1199]coordinates: Latitude Longitude\nlong_name: Band 31 emissivity\nEmis_31: Emis_31 data * scale_factor + add_offset\nscale_factor_err: 0.0000000000000000\nadd_offset_err: 0.0000000000000000\ncalibrated_nt: 5\norig_scale_factor: 0.0020000000000000000\norig_add_offset: 0.48999999999999999\n_FillValue: 0\nvalid_min: 0.4920000136\nvalid_max: 1.000000000\nNumber_Type_Orig: uint8\ngrid_mapping: MODIS_Grid_Daily_1km_LST_eos_cf_projection
#> 17                                                     Array of 32 bit Reals [time = 0..6844][YDim = 0..1199][XDim = 0..1199]coordinates: Latitude Longitude\nlong_name: Band 32 emissivity\nEmis_32: Emis_32 data * scale_factor + add_offset\nscale_factor_err: 0.0000000000000000\nadd_offset_err: 0.0000000000000000\ncalibrated_nt: 5\norig_scale_factor: 0.0020000000000000000\norig_add_offset: 0.48999999999999999\n_FillValue: 0\nvalid_min: 0.4920000136\nvalid_max: 1.000000000\nNumber_Type_Orig: uint8\ngrid_mapping: MODIS_Grid_Daily_1km_LST_eos_cf_projection
#> 18                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   Array of 64 bit Reals [time = 0..6844]units: days since 2000-01-01 00:00
#>     extractable_w_opendapr
#> 13             extractable
#> 14             extractable
#> 15             extractable
#> 16             extractable
#> 17             extractable
#> 18 automatically extracted

In odr_get_url(), the parameter variables enables to restrict the data to download to only specific variables. By default it is set to NULL, which means that all the available variables for the specified collection are downloaded. Specifying variables will make the data to download lighter.

End Note

In our example, we specify some variables for each collection.

## Get the URLs of MODIS Terra LST daily
urls_mod11a1 <- odr_get_url(
  collection = "MOD11A1.006",
  variables = c("LST_Day_1km","LST_Night_1km","QC_Day","QC_Night"),  # get the variables available with : odr_list_variables("MOD11A1.006") ; or set to NULL (defaults) to download all the variables
  roi = roi,
  time_range = time_range
 )
#> Building the URLs...
#> Note : messages above ('although coordinates are longitude/latitude, st_intersection assumes that they are planar' and 'attribute variables are assumed to be spatially constant throughout all geometries') are not errors
#> OK

## Get the URLs of GPM daily
urls_gpm <- odr_get_url(
  collection = "GPM_3IMERGDF.06",
  variables = c("precipitationCal","precipitationCal_cnt"),  # get the variables available with : odr_list_variables("GPM_L3/GPM_3IMERGDF.06")
  roi = roi,
  time_range = time_range
 )
#> Building the URLs...
#> OK

## Get the URLs of SMAP 3-days
urls_smap <- odr_get_url(
  collection = "SPL3SMP_E.003",
  variables = c("Soil_Moisture_Retrieval_Data_AM_soil_moisture","Soil_Moisture_Retrieval_Data_AM_retrieval_qual_flag","Soil_Moisture_Retrieval_Data_PM_soil_moisture_pm","Soil_Moisture_Retrieval_Data_PM_retrieval_qual_flag_pm"),     # get the variables available with : odr_list_variables("SMAP/SPL3SMP_E.003")
  roi = roi,
  time_range = time_range
 )
#> Building the URLs...
#> OK


nrow(urls_mod11a1)
#> [1] 1
head(urls_mod11a1,3)
#>   time_start                               name
#> 1 2017-01-01 MOD11A1.006.2017001_2017030.h17v08
#>                                                                                                                                                                                                                                                                                                                     url
#> 1 https://opendap.cr.usgs.gov/opendap/hyrax/MOD11A1.006/h17v08.ncml.nc4?MODIS_Grid_Daily_1km_LST_eos_cf_projection,LST_Day_1km[6093:6122][55:140][512:560],LST_Night_1km[6093:6122][55:140][512:560],QC_Day[6093:6122][55:140][512:560],QC_Night[6093:6122][55:140][512:560],time[6093:6122],YDim[55:140],XDim[512:560]
#>                                             destfile
#> 1 MOD11A1.006/MOD11A1.006.2017001_2017030.h17v08.nc4

nrow(urls_gpm)
#> [1] 30
head(urls_gpm,3)
#>   time_start                                              name
#> 1 2017-01-01 3B-DAY.MS.MRG.3IMERG.20170101-S000000-E235959.V06
#> 2 2017-01-02 3B-DAY.MS.MRG.3IMERG.20170102-S000000-E235959.V06
#> 3 2017-01-03 3B-DAY.MS.MRG.3IMERG.20170103-S000000-E235959.V06
#>                                                                                                                                                                                                                                                                   url
#> 1 https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDF.06/2017/01/3B-DAY.MS.MRG.3IMERG.20170101-S000000-E235959.V06.nc4.nc4?precipitationCal[0:0][1742:1746][989:996],precipitationCal_cnt[0:0][1742:1746][989:996],time[0:0],lon[1742:1746],lat[989:996]
#> 2 https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDF.06/2017/01/3B-DAY.MS.MRG.3IMERG.20170102-S000000-E235959.V06.nc4.nc4?precipitationCal[0:0][1742:1746][989:996],precipitationCal_cnt[0:0][1742:1746][989:996],time[0:0],lon[1742:1746],lat[989:996]
#> 3 https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDF.06/2017/01/3B-DAY.MS.MRG.3IMERG.20170103-S000000-E235959.V06.nc4.nc4?precipitationCal[0:0][1742:1746][989:996],precipitationCal_cnt[0:0][1742:1746][989:996],time[0:0],lon[1742:1746],lat[989:996]
#>                                                                destfile
#> 1 GPM_3IMERGDF.06/3B-DAY.MS.MRG.3IMERG.20170101-S000000-E235959.V06.nc4
#> 2 GPM_3IMERGDF.06/3B-DAY.MS.MRG.3IMERG.20170102-S000000-E235959.V06.nc4
#> 3 GPM_3IMERGDF.06/3B-DAY.MS.MRG.3IMERG.20170103-S000000-E235959.V06.nc4

nrow(urls_smap)
#> [1] 30
head(urls_smap,3)
#>   time_start                               name
#> 1 2017-01-01 SMAP_L3_SM_P_E_20170101_R16510_001
#> 2 2017-01-02 SMAP_L3_SM_P_E_20170102_R16510_001
#> 3 2017-01-03 SMAP_L3_SM_P_E_20170103_R16510_001
#>                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                    url
#> 1 https://n5eil02u.ecs.nsidc.org/opendap/SMAP/SPL3SMP_E.003/2017.01.01/SMAP_L3_SM_P_E_20170101_R16510_001.h5.nc4?Soil_Moisture_Retrieval_Data_AM_soil_moisture[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_retrieval_qual_flag[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_soil_moisture_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_retrieval_qual_flag_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_longitude[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_latitude[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_longitude_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_latitude_pm[678:1:688][1866:1:1871]
#> 2 https://n5eil02u.ecs.nsidc.org/opendap/SMAP/SPL3SMP_E.003/2017.01.02/SMAP_L3_SM_P_E_20170102_R16510_001.h5.nc4?Soil_Moisture_Retrieval_Data_AM_soil_moisture[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_retrieval_qual_flag[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_soil_moisture_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_retrieval_qual_flag_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_longitude[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_latitude[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_longitude_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_latitude_pm[678:1:688][1866:1:1871]
#> 3 https://n5eil02u.ecs.nsidc.org/opendap/SMAP/SPL3SMP_E.003/2017.01.03/SMAP_L3_SM_P_E_20170103_R16510_001.h5.nc4?Soil_Moisture_Retrieval_Data_AM_soil_moisture[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_retrieval_qual_flag[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_soil_moisture_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_retrieval_qual_flag_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_longitude[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_AM_latitude[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_longitude_pm[678:1:688][1866:1:1871],Soil_Moisture_Retrieval_Data_PM_latitude_pm[678:1:688][1866:1:1871]
#>                                               destfile
#> 1 SPL3SMP_E.003/SMAP_L3_SM_P_E_20170101_R16510_001.nc4
#> 2 SPL3SMP_E.003/SMAP_L3_SM_P_E_20170102_R16510_001.nc4
#> 3 SPL3SMP_E.003/SMAP_L3_SM_P_E_20170103_R16510_001.nc4

Download the data

Now we download the data with the function odr_download_data().

Destination file for each dataset is specified in the column destfile of the dataframes urls_mod11a1, urls_gpm and urls_smap. The destination file is specified by default but it can of course be modified.

Setting the argument parallel to TRUE will parallelize - therefore fasten - the download in case their are numerous datasets to download

df_to_dl <- rbind(urls_mod11a1,urls_gpm,urls_smap)
res_dl <- odr_download_data(df_to_dl,source="earthdata",parallel = TRUE)
#> 61  datasets in total :  61  already downloaded and  0  datasets to download
#> OK

print(str(res_dl))
#> 'data.frame':    61 obs. of  7 variables:
#>  $ time_start: Date, format: "2017-01-01" "2017-01-01" ...
#>  $ name      : chr  "MOD11A1.006.2017001_2017030.h17v08" "3B-DAY.MS.MRG.3IMERG.20170101-S000000-E235959.V06" "3B-DAY.MS.MRG.3IMERG.20170102-S000000-E235959.V06" "3B-DAY.MS.MRG.3IMERG.20170103-S000000-E235959.V06" ...
#>  $ url       : chr  "https://opendap.cr.usgs.gov/opendap/hyrax/MOD11A1.006/h17v08.ncml.nc4?MODIS_Grid_Daily_1km_LST_eos_cf_projectio"| __truncated__ "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDF.06/2017/01/3B-DAY.MS.MRG.3IMERG.20170101-S0000"| __truncated__ "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDF.06/2017/01/3B-DAY.MS.MRG.3IMERG.20170102-S0000"| __truncated__ "https://gpm1.gesdisc.eosdis.nasa.gov/opendap/GPM_L3/GPM_3IMERGDF.06/2017/01/3B-DAY.MS.MRG.3IMERG.20170103-S0000"| __truncated__ ...
#>  $ destfile  : chr  "MOD11A1.006/MOD11A1.006.2017001_2017030.h17v08.nc4" "GPM_3IMERGDF.06/3B-DAY.MS.MRG.3IMERG.20170101-S000000-E235959.V06.nc4" "GPM_3IMERGDF.06/3B-DAY.MS.MRG.3IMERG.20170102-S000000-E235959.V06.nc4" "GPM_3IMERGDF.06/3B-DAY.MS.MRG.3IMERG.20170103-S000000-E235959.V06.nc4" ...
#>  $ fileDl    : logi  TRUE TRUE TRUE TRUE TRUE TRUE ...
#>  $ fileSize  : num  287343 34510 34516 34518 34510 ...
#>  $ dlStatus  : num  3 3 3 3 3 3 3 3 3 3 ...
#> NULL

Import the data in R

Importing the data in R goes beyond the scope of this package. However, here we provide some tips and tricks to do the job for further data analysis in R !

Various packages and related classes can be used to read the datacubes downloaded through OPeNDAP (in NetCDF, ascii or json format). If raster is surely the most famous class for raster objects, many packages facilitate the use of spatiotemporal data cubes in formats such as those proposed through opendapr (e.g. NetCDF). For instance, MODIS or VIIRS products can be imported as a stars object from the excellent stars package for data cubes manipulation. All the data can also be imported as ncdf4 objects using e.g. the ncdf4 package, or RasterLayer of the raster package.

Important note regarding the import of the data in R

In any case, care must be taken when importing data that was downloaded through the OPeNDAP data providers servers. Depending on the collection, some “issues” were raised. These issues are independant from opendapr : they result most of time of a kind of lack of full implementation of the OPeNDAP framework by the data providers. These issues are :

  • for MODIS and VIIRS collections : CRS has to be provided
  • for GPM collections : CRS has to be provided + data have to be flipped
  • for SMAP collections : CRS + bounding coordinates of the data have to be provided

These issues can easily be dealt at the import phase in R. Below we propose some functions that include the processings that have to be done at the data import phase in order to open the data as raster objects. (argument destfiles is the path to a dataset downloaded with opendapr - output of odr_get_url()$destfile - and variable is the name of a variable to import).

Functions to import the data in R

Import MODIS / VIIRS

As a RasterLayer object :

## Functions to import MODIS and VIIRS products as RasterLayer object
# In case the ROI covers one single MODIS tile :
require(raster)
.import_modis_onetile <- function(destfiles,variable){
  rasts <- destfiles %>%
    raster::brick(.,varname=variable,crs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")
  return(rasts)
}

# In case the ROI covers multiple MODIS tiles :
.import_modis_moretiles <- function(destfiles,variable){
  rasts <- destfiles %>%
    purrr::map(~raster::brick(.,varname=variable,crs="+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs")) %>%
    do.call(merge,.)
  return(rasts)
}

As a stars object :

We can also import the same MODIS time series as a stars object. Here the interesting point is that all the dimensions are imported at once.

Import GPM

As a RasterLayer object :

## Function to import GPM products as RasterLayer object
require(purrr)
.import_gpm <- function(destfiles,variable){
  rasts <- destfiles %>%
    purrr::map(~raster(., varname = variable,crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ")) %>%
    raster::brick() %>%
    raster::t() %>%
    raster::flip("y") %>%
    raster::flip("x")
  return(rasts)
}

Import SMAP

As a RasterLayer object :

We set-up the missing bounding coordinates. For this we use the function odr_get_opt_param().

## Function to import SMAP products as RasterLayer object
require(purrr)
require(ncdf4)
require(raster)

smap_sp_bound <- opendapr::odr_get_opt_param(roi = roi, collection = "SPL3SMP_E.003")$roiSpatialBound$`1`

.import_smap <- function(destfiles,variable,smap_sp_bound){
 rasts <- destfiles %>%
   purrr::map(~ncdf4::nc_open(.)) %>%
   purrr::map(~ncdf4::ncvar_get(., "Soil_Moisture_Retrieval_Data_AM_soil_moisture")) %>%
   purrr::map(~raster(t(.), ymn=smap_sp_bound[1], ymx=smap_sp_bound[2], xmn=smap_sp_bound[3], xmx=smap_sp_bound[4], crs="+proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0")) %>%  # EPSG : 6933
   raster::brick()
  return(rasts)
}

Import the data from our example

So let’s import the data from our example, using the functions provided below :

# import MOD11A1.006
mod11a1_rast_day <- .import_modis_onetile(urls_mod11a1$destfile, "LST_Day_1km")
# import GPM
gpm_rast_precipitationcal <- .import_gpm(urls_gpm$destfile, "precipitationCal")
# import SPL3SMP_E.003
smap_rast_sm_am <- .import_smap(urls_smap$destfile, "Soil_Moisture_Retrieval_Data_AM_soil_moisture",smap_sp_bound)

mod11a1_rast_day
#> class       : RasterBrick 
#> dimensions  : 86, 48, 4128, 30  (nrow, ncol, ncell, nlayers)
#> resolution  : 926.6254, 926.6254  (x, y)
#> extent      : -637981.6, -593503.6, 981759.6, 1061449  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs 
#> data source : /home/ptaconet/opendapr/vignettes/MOD11A1.006/MOD11A1.006.2017001_2017030.h17v08.nc4 
#> names       : X2017.01.01, X2017.01.02, X2017.01.03, X2017.01.04, X2017.01.05, X2017.01.06, X2017.01.07, X2017.01.08, X2017.01.09, X2017.01.10, X2017.01.11, X2017.01.12, X2017.01.13, X2017.01.14, X2017.01.15, ... 
#> Date        : 2017-01-01, 2017-01-30 (min, max)
#> varname     : LST_Day_1km
gpm_rast_precipitationcal
#> class       : RasterBrick 
#> dimensions  : 8, 5, 40, 30  (nrow, ncol, ncell, nlayers)
#> resolution  : 0.1000023, 0.09999956  (x, y)
#> extent      : -5.800004, -5.299993, 8.900002, 9.699998  (xmin, xmax, ymin, ymax)
#> coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
#> data source : in memory
#> names       :     layer.1,     layer.2,     layer.3,     layer.4,     layer.5,     layer.6,     layer.7,     layer.8,     layer.9,    layer.10,    layer.11,    layer.12,    layer.13,    layer.14,    layer.15, ... 
#> min values  : 0.000000000, 0.000000000, 0.000000000, 0.000000000, 0.000000000, 0.000000000, 0.000000000, 0.000000000, 0.000000000, 0.015777921, 0.000000000, 0.021967307, 0.002506131, 0.000000000, 0.000000000, ... 
#> max values  :  0.00000000,  1.00454104,  0.01239155,  0.00000000,  0.00000000,  0.00000000,  0.00000000,  0.02110096,  0.34920487,  0.42241478,  2.99069405,  5.54364347,  0.19833815,  0.21644105,  1.17732072, ...
smap_rast_sm_am
#> class       : RasterBrick 
#> dimensions  : 11, 6, 66, 30  (nrow, ncol, ncell, nlayers)
#> resolution  : 7506.713, 8189.135  (x, y)
#> extent      : -563003.5, -517963.2, 1121503, 1211583  (xmin, xmax, ymin, ymax)
#> coord. ref. : +proj=cea +lon_0=0 +lat_ts=30 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0 
#> data source : in memory
#> names       :    layer.1,    layer.2,    layer.3,    layer.4,    layer.5,    layer.6,    layer.7,    layer.8,    layer.9,   layer.10,   layer.11,   layer.12,   layer.13,   layer.14,   layer.15, ... 
#> min values  :         NA, 0.09996461,         NA,         NA, 0.09471568,         NA, 0.05590146,         NA,         NA, 0.09872495,         NA,         NA, 0.09351711,         NA, 0.06066480, ... 
#> max values  :         NA,  0.1452434,         NA,         NA,  0.1560785,         NA,  0.1505869,         NA,         NA,  0.1704942,         NA,         NA,  0.1503100,         NA,  0.1584005, ...

Plot the data

Let’s finally plot the data !

(Note that only the first 16 dates are plotted here-under)

# Land surface temperature
mod11a1_rast_day <- projectRaster(mod11a1_rast_day,crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ")
plot(mod11a1_rast_day)

# Precipitation
names(gpm_rast_precipitationcal) <- urls_gpm$time_start
plot(gpm_rast_precipitationcal)

# Soil moisture
smap_rast_sm_am <- projectRaster(smap_rast_sm_am,crs = "+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 ")
names(smap_rast_sm_am) <- urls_smap$time_start
plot(smap_rast_sm_am)