Geospatial Datacube Processing with R
This notebook gives a short introduction to using R with OGC WCS and WCPS, based on examples hosted in a rasdaman datacube service. It is based on a Bachelor Thesis done by Alina Amanbayeva at Constructor University done with the Large-Scale Scientific Information Systems Research Group.
How to integrate WCS requests into an R terminal¶
You can install required R libraries using install.packages("name_of_package_listed_in_CRAN") function.
List of all packages can be found here: https://cran.r-project.org/web/packages/available_packages_by_name.html
To save time I created a package already containing all the packages needed to download. In order to use it run:
In [3]:
#load pre-installed packages needed
install.packages("devtools")
install.packages("remotes")
library(devtools)
remotes::install_github("aamanbayev/CoverageProcessingR", dependencies = TRUE)
devtools::install_github("ARPASMR/myCubeR@HEAD")
library(CoverageProcessingR)
Creation of an interface in R¶
In order to operate on a Web Coverage Service, you need to create an interface in R to this WCS. This is done with the class WCSClient, as follows:
In [4]:
wcs_client <- ows4R::WCSClient$new("https://ows.rasdaman.org/rasdaman/ows", "2.1.0", logger = "INFO")
[ows4R][INFO] OWSGetCapabilities - Fetching https://ows.rasdaman.org/rasdaman/ows?service=WCS&version=2.1.0&request=GetCapabilities |======================================================================| 100%
Checking the connection¶
In [5]:
caps <- tryCatch(
wcs_client$getCapabilities(),
error = function(e) {
message("Error: Could not connect to WCS server.")
return(NULL)
}
) #Parse a message in case if it's successful
if (!is.null(caps)) {
message("Connection to WCS server successful.")
# Do something with the capabilities
} else {
message("Connection to WCS server failed.")
}#Let the user know if connection failed
Connection to WCS server successful.
Disclosing all available coverages¶
Print available coverage IDs taken from the result of GetCoverageSummaries() function
In [7]:
coverage_summaries <- caps$getCoverageSummaries()
coverage_ids <- sapply(coverage_summaries, function(summary) summary$CoverageId)
print(coverage_ids)
[1] "AverageChloroColor" [2] "AverageChloroColorScaled" [3] "AvgLandTemp" [4] "AvgTemperatureColor" [5] "AvgTemperatureColorScaled" [6] "Bavaria_50_DSM" [7] "BlueMarbleCov" [8] "Germany_DTM" [9] "NIR" [10] "NN3_1" [11] "NN3_2" [12] "NN3_3" [13] "NN3_4" [14] "RadianceColorScaled" [15] "S2_L2A_32631_B01_60m" [16] "S2_L2A_32631_B03_10m" [17] "S2_L2A_32631_B04_10m" [18] "S2_L2A_32631_B08_10m" [19] "S2_L2A_32631_B12_20m" [20] "S2_L2A_32631_TCI_60m" [21] "S2_federation_demo" [22] "Temperature4D" [23] "climate_cloud" [24] "climate_earth" [25] "land_cover_class__esa_test_6" [26] "mean_summer_airtemp" [27] "multiband" [28] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226" [29] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_128" [30] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_16" [31] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_2" [32] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_256" [33] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_32" [34] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_4" [35] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_64" [36] "nchc.earthserver.xyz:7000:AER_DLA_DEM_4326_2000_DMC_20010814_20061226_8" [37] "nchc.earthserver.xyz:7000:AverageChlorophyll_new" [38] "nchc.earthserver.xyz:7000:AverageTemperature" [39] "nchc.earthserver.xyz:7000:LANDUSE_50Q_20200101_20210101" [40] "nchc.earthserver.xyz:7000:LANDUSE_50Q_20200101_20210101_2" [41] "nchc.earthserver.xyz:7000:LANDUSE_50Q_20200101_20210101_4" [42] "nchc.earthserver.xyz:7000:LANDUSE_50Q_20200101_20210101_8" [43] "nchc.earthserver.xyz:7000:LANDUSE_51Q_20200101_20210101" [44] "nchc.earthserver.xyz:7000:LANDUSE_51Q_20200101_20210101_2" [45] "nchc.earthserver.xyz:7000:LANDUSE_51Q_20200101_20210101_4" [46] "nchc.earthserver.xyz:7000:LANDUSE_51Q_20200101_20210101_8" [47] "nchc.earthserver.xyz:7000:LANDUSE_51R_20200101_20210101" [48] "nchc.earthserver.xyz:7000:LANDUSE_51R_20200101_20210101_2" [49] "nchc.earthserver.xyz:7000:LANDUSE_51R_20200101_20210101_4" [50] "nchc.earthserver.xyz:7000:LANDUSE_51R_20200101_20210101_8" [51] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101" [52] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_128" [53] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_16" [54] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_2" [55] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_256" [56] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_32" [57] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_4" [58] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_64" [59] "nchc.earthserver.xyz:7000:LIDAR_DLA_DEM_4326_2000_Leica_20110101_20161101_8" [60] "nchc.earthserver.xyz:7000:MSMlvltest" [61] "nchc.earthserver.xyz:7000:MSMlvltest_2" [62] "nchc.earthserver.xyz:7000:MSMlvltest_4" [63] "nchc.earthserver.xyz:7000:TW_DLA_20010814_20061226_20M_3826_DEM" [64] "nchc.earthserver.xyz:7000:TW_DLA_20100101_20191101_20M_3826_DEM" [65] "nchc.earthserver.xyz:7000:TW_DLA_20110101_20161101_20M_3826_DEM" [66] "nchc.earthserver.xyz:7000:m_msm_lvl_2021" [67] "nchc.earthserver.xyz:7000:m_msm_lvl_2021_2" [68] "nchc.earthserver.xyz:7000:m_msm_lvl_2021_4" [69] "nchc.earthserver.xyz:7000:m_msm_sfc_2021" [70] "nchc.earthserver.xyz:7000:m_msm_sfc_2021_2" [71] "nchc.earthserver.xyz:7000:m_msm_sfc_2021_4" [72] "nchc.earthserver.xyz:7000:mean_summer_airtemp" [73] "nchc.earthserver.xyz:7000:sample" [74] "nchc.earthserver.xyz:7000:sfc_allp" [75] "nchc.earthserver.xyz:7000:sfc_allp_2" [76] "nchc.earthserver.xyz:7000:sfc_allp_4" [77] "nchc.earthserver.xyz:7000:sfc_test1" [78] "nchc.earthserver.xyz:7000:sfc_test1_2" [79] "nchc.earthserver.xyz:7000:sfc_test1_4" [80] "nchc.earthserver.xyz:7000:sfc_test_4" [81] "nchc.earthserver.xyz:7000:taiwan_air_quality_monitoring_data" [82] "nchc.earthserver.xyz:7000:taiwan_air_quality_monitoring_data_benzene" [83] "nchc.earthserver.xyz:7000:taiwan_air_quality_monitoring_data_co" [84] "nchc.earthserver.xyz:7000:taiwan_air_quality_monitoring_data_no2" [85] "nchc.earthserver.xyz:7000:taiwan_air_quality_monitoring_data_o3" [86] "nchc.earthserver.xyz:7000:taiwan_air_quality_monitoring_data_pm2dot5" [87] "nchc.earthserver.xyz:7000:taiwan_air_quality_monitoring_data_so2"
In [8]:
chla <- caps$findCoverageSummaryById("AverageChloroColorScaled", exact = T)
print(chla)
<WCSCoverageSummary>
Inherits from: <OGCAbstractObject>
Public:
attrs: list
BoundingBox: list
clone: function (deep = FALSE)
CoverageId: AverageChloroColorScaled
CoverageSubtype: ReferenceableGridCoverage
CoverageSubtypeParent: NULL
defaults: list
element: AbstractObject
encode: function (addNS = TRUE, geometa_validate = TRUE, geometa_inspire = FALSE,
ERROR: function (text)
getBoundingBox: function ()
getClass: function ()
getClassName: function ()
getCoverage: function (bbox = NULL, crs = NULL, time = NULL, elevation = NULL,
getCoverageStack: function (time = NULL, elevation = NULL, bbox = NULL, filename_handler = NULL,
getDescription: function ()
getDimensions: function ()
getId: function ()
getNamespaceDefinition: function (recursive = FALSE)
getSubtype: function ()
getSubtypeParent: function ()
getWGS84BoundingBox: function ()
INFO: function (text)
initialize: function (xmlObj, capabilities, serviceVersion, owsVersion, logger = NULL)
isFieldInheritedFrom: function (field)
logger: function (type, text)
loggerType: INFO
namespace: OWSNamespace, R6
verbose.debug: FALSE
verbose.info: TRUE
WARN: function (text)
WGS84BoundingBox: list
wrap: FALSE
Private:
capabilities: WCSCapabilities, OWSCapabilities, OGCAbstractObject, R6
description: NULL
dimensions: NULL
fetchCoverageSummary: function (xmlObj, serviceVersion, owsVersion)
fromComplexTypes: function (value)
owsVersion: 2.0
system_fields: verbose.info verbose.debug loggerType wrap element names ...
url: https://ows.rasdaman.org/rasdaman/ows
version: 2.1.0
xmlElement: AbstractObject
xmlExtraNamespaces: NULL
xmlNamespacePrefix: OWS
xmlNodeToCharacter: function (x, ..., indent = "", tagSeparator = "\n")
Get coverage axis labels with order according to coverage's Coordinate Reference System (CRS)¶
In [10]:
chla_dims <- chla$getDimensions()
for (dim in chla_dims) {
print(dim$label)
[ows4R][INFO] WCSDescribeCoverage - Fetching https://ows.rasdaman.org/rasdaman/ows?service=WCS&version=2.1.0&coverageId=AverageChloroColorScaled&request=DescribeCoverage |======================================================================| 100% [ows4R][INFO] WCSCoverageSummary - Fetching Coverage envelope dimensions by CRS interpretation [ows4R][INFO] WCSCoverageSummary - Try to parse CRS from 'http://www.opengis.net/def/crs/OGC/0/AnsiDate' [ows4R][INFO] WCSCoverageSummary - Try to parse CRS from 'http://www.opengis.net/def/crs/EPSG/0/4326' [1] "ansi" [1] "Lat" [1] "Lon"
Displaying the result of the describeCoverage function¶
which includes grid's upper and lower limits,and offsetvectors in case you would like the result to be more user friendly getDescription() function can be used
In [14]:
chla_des <- wcs_client$describeCoverage("AverageChloroColorScaled")
print(chla_des)
[ows4R][INFO] WCSClient - Fetching coverageSummary description for 'AverageChloroColorScaled' ... <WCSCoverageDescription> ....|-- boundedBy [srsName=http://crs.rasdaman.com/def/crs-compound?1=http://crs.rasdaman.com/def/crs/OGC/0/AnsiDate&2=http://crs.rasdaman.com/def/crs/EPSG/0/4326,axisLabels=ansi Lat Lon,uomLabels=d degree degree,srsDimension=3] <GMLEnvelope> ........|-- lowerCorner: "2002-07-01T00:00:00.000Z" -90 -180 ........|-- upperCorner: "2015-05-01T00:00:00.000Z" 90 180 ....|-- domainSet [dimension=3,gmlrgrid:id=AverageChloroColorScaled-grid] <GMLReferenceableGridByVectors> ........|-- limits <GMLGridEnvelope> ............|-- low ................|-- value: 0 0 0 ............|-- high ................|-- value: 154 449 899 ........|-- axisLabels ............|-- value: ansi Lat Lon ........|-- origin [gml:id=AverageChloroColorScaled-point,srsName=http://crs.rasdaman.com/def/crs-compound?1=http://crs.rasdaman.com/def/crs/OGC/0/AnsiDate&2=http://crs.rasdaman.com/def/crs/EPSG/0/4326] <GMLPoint> ............|-- pos: "2002-07-01T00:00:00.000Z" 89.8 -179.8 ........|-- generalGridAxis <GMLGeneralGridAxis> ............|-- offsetVector: 1 0 0 ............|-- coefficients: 2002-07-01T00:00:00.000Z 2002-08-01T00:00:00.000Z 2002-09-01T00:00:00.000Z 2002-10-01T00:00:00.000Z 2002-11-01T00:00:00.000Z 2002-12-01T00:00:00.000Z 2003-01-01T00:00:00.000Z 2003-02-01T00:00:00.000Z 2003-03-01T00:00:00.000Z 2003-04-01T00:00:00.000Z 2003-05-01T00:00:00.000Z 2003-06-01T00:00:00.000Z 2003-07-01T00:00:00.000Z 2003-08-01T00:00:00.000Z 2003-09-01T00:00:00.000Z 2003-10-01T00:00:00.000Z 2003-11-01T00:00:00.000Z 2003-12-01T00:00:00.000Z 2004-01-01T00:00:00.000Z 2004-02-01T00:00:00.000Z 2004-03-01T00:00:00.000Z 2004-04-01T00:00:00.000Z 2004-05-01T00:00:00.000Z 2004-06-01T00:00:00.000Z 2004-07-01T00:00:00.000Z 2004-08-01T00:00:00.000Z 2004-09-01T00:00:00.000Z 2004-10-01T00:00:00.000Z 2004-11-01T00:00:00.000Z 2004-12-01T00:00:00.000Z 2005-01-01T00:00:00.000Z 2005-02-01T00:00:00.000Z 2005-03-01T00:00:00.000Z 2005-04-01T00:00:00.000Z 2005-05-01T00:00:00.000Z 2005-06-01T00:00:00.000Z 2005-07-01T00:00:00.000Z 2005-08-01T00:00:00.000Z 2005-09-01T00:00:00.000Z 2005-10-01T00:00:00.000Z 2005-11-01T00:00:00.000Z 2005-12-01T00:00:00.000Z 2006-01-01T00:00:00.000Z 2006-02-01T00:00:00.000Z 2006-03-01T00:00:00.000Z 2006-04-01T00:00:00.000Z 2006-05-01T00:00:00.000Z 2006-06-01T00:00:00.000Z 2006-07-01T00:00:00.000Z 2006-08-01T00:00:00.000Z 2006-09-01T00:00:00.000Z 2006-10-01T00:00:00.000Z 2006-11-01T00:00:00.000Z 2006-12-01T00:00:00.000Z 2007-01-01T00:00:00.000Z 2007-02-01T00:00:00.000Z 2007-03-01T00:00:00.000Z 2007-04-01T00:00:00.000Z 2007-05-01T00:00:00.000Z 2007-06-01T00:00:00.000Z 2007-07-01T00:00:00.000Z 2007-08-01T00:00:00.000Z 2007-09-01T00:00:00.000Z 2007-10-01T00:00:00.000Z 2007-11-01T00:00:00.000Z 2007-12-01T00:00:00.000Z 2008-01-01T00:00:00.000Z 2008-02-01T00:00:00.000Z 2008-03-01T00:00:00.000Z 2008-04-01T00:00:00.000Z 2008-05-01T00:00:00.000Z 2008-06-01T00:00:00.000Z 2008-07-01T00:00:00.000Z 2008-08-01T00:00:00.000Z 2008-09-01T00:00:00.000Z 2008-10-01T00:00:00.000Z 2008-11-01T00:00:00.000Z 2008-12-01T00:00:00.000Z 2009-01-01T00:00:00.000Z 2009-02-01T00:00:00.000Z 2009-03-01T00:00:00.000Z 2009-04-01T00:00:00.000Z 2009-05-01T00:00:00.000Z 2009-06-01T00:00:00.000Z 2009-07-01T00:00:00.000Z 2009-08-01T00:00:00.000Z 2009-09-01T00:00:00.000Z 2009-10-01T00:00:00.000Z 2009-11-01T00:00:00.000Z 2009-12-01T00:00:00.000Z 2010-01-01T00:00:00.000Z 2010-02-01T00:00:00.000Z 2010-03-01T00:00:00.000Z 2010-04-01T00:00:00.000Z 2010-05-01T00:00:00.000Z 2010-06-01T00:00:00.000Z 2010-07-01T00:00:00.000Z 2010-08-01T00:00:00.000Z 2010-09-01T00:00:00.000Z 2010-10-01T00:00:00.000Z 2010-11-01T00:00:00.000Z 2010-12-01T00:00:00.000Z 2011-01-01T00:00:00.000Z 2011-02-01T00:00:00.000Z 2011-03-01T00:00:00.000Z 2011-04-01T00:00:00.000Z 2011-05-01T00:00:00.000Z 2011-06-01T00:00:00.000Z 2011-07-01T00:00:00.000Z 2011-08-01T00:00:00.000Z 2011-09-01T00:00:00.000Z 2011-10-01T00:00:00.000Z 2011-11-01T00:00:00.000Z 2011-12-01T00:00:00.000Z 2012-01-01T00:00:00.000Z 2012-02-01T00:00:00.000Z 2012-03-01T00:00:00.000Z 2012-04-01T00:00:00.000Z 2012-05-01T00:00:00.000Z 2012-06-01T00:00:00.000Z 2012-07-01T00:00:00.000Z 2012-08-01T00:00:00.000Z 2012-09-01T00:00:00.000Z 2012-10-01T00:00:00.000Z 2012-11-01T00:00:00.000Z 2012-12-01T00:00:00.000Z 2013-01-01T00:00:00.000Z 2013-02-01T00:00:00.000Z 2013-03-01T00:00:00.000Z 2013-04-01T00:00:00.000Z 2013-05-01T00:00:00.000Z 2013-06-01T00:00:00.000Z 2013-07-01T00:00:00.000Z 2013-08-01T00:00:00.000Z 2013-09-01T00:00:00.000Z 2013-10-01T00:00:00.000Z 2013-11-01T00:00:00.000Z 2013-12-01T00:00:00.000Z 2014-01-01T00:00:00.000Z 2014-02-01T00:00:00.000Z 2014-03-01T00:00:00.000Z 2014-04-01T00:00:00.000Z 2014-05-01T00:00:00.000Z 2014-06-01T00:00:00.000Z 2014-07-01T00:00:00.000Z 2014-08-01T00:00:00.000Z 2014-09-01T00:00:00.000Z 2014-10-01T00:00:00.000Z 2014-11-01T00:00:00.000Z 2014-12-01T00:00:00.000Z 2015-01-01T00:00:00.000Z 2015-02-01T00:00:00.000Z 2015-03-01T00:00:00.000Z 2015-04-01T00:00:00.000Z 2015-05-01T00:00:00.000Z ............|-- gridAxesSpanned ................|-- value: ansi ............|-- sequenceRule [axisOrder=+1] ................|-- value: Linear ........|-- generalGridAxis <GMLGeneralGridAxis> ............|-- offsetVector: 0 -0.4 0 ............|-- coefficients ............|-- gridAxesSpanned ................|-- value: Lat ............|-- sequenceRule [axisOrder=+1] ................|-- value: Linear ........|-- generalGridAxis <GMLGeneralGridAxis> ............|-- offsetVector: 0 0 0.4 ............|-- coefficients ............|-- gridAxesSpanned ................|-- value: Lon ............|-- sequenceRule [axisOrder=+1] ................|-- value: Linear ....|-- coverageFunction <GMLGridFunction> ........|-- sequenceRule [axisOrder=+2 +3 +1] ............|-- value: Linear ........|-- startPoint ............|-- value: 0 0 0 ....|-- rangeType <SWEDataRecord> ........|-- field [definition=http://www.opengis.net/def/dataType/OGC/0/unsignedByte]<SWEQuantity> ............|-- label ................|-- value: Red ............|-- description ............|-- uom [code=10^0] ............|-- constraint ........|-- field [definition=http://www.opengis.net/def/dataType/OGC/0/unsignedByte]<SWEQuantity> ............|-- label ................|-- value: Green ............|-- description ............|-- uom [code=10^0] ............|-- constraint ........|-- field [definition=http://www.opengis.net/def/dataType/OGC/0/unsignedByte]<SWEQuantity> ............|-- label ................|-- value: Blue ............|-- description ............|-- uom [code=10^0] ............|-- constraint ....|-- metadata ....|-- CoverageId ........|-- value: AverageChloroColorScaled ....|-- ServiceParameters <ISOElementSequence> ........|-- CoverageSubtype: ReferenceableGridCoverage ........|-- nativeFormat: application/octet-stream
Get the list of coefficients for irregular axis (time is an irregular axis of this coverage)¶
In [8]:
for (dim in chla_dims) {
if (dim$type == "temporal") {
print(dim$coefficients)
}
}
[,1] [1,] "2002-07-01T00:00:00.000Z" [2,] "2002-08-01T00:00:00.000Z" [3,] "2002-09-01T00:00:00.000Z" [4,] "2002-10-01T00:00:00.000Z" [5,] "2002-11-01T00:00:00.000Z" [6,] "2002-12-01T00:00:00.000Z" [7,] "2003-01-01T00:00:00.000Z" [8,] "2003-02-01T00:00:00.000Z" [9,] "2003-03-01T00:00:00.000Z" [10,] "2003-04-01T00:00:00.000Z" [11,] "2003-05-01T00:00:00.000Z" [12,] "2003-06-01T00:00:00.000Z" [13,] "2003-07-01T00:00:00.000Z" [14,] "2003-08-01T00:00:00.000Z" [15,] "2003-09-01T00:00:00.000Z" [16,] "2003-10-01T00:00:00.000Z" [17,] "2003-11-01T00:00:00.000Z" [18,] "2003-12-01T00:00:00.000Z" [19,] "2004-01-01T00:00:00.000Z" [20,] "2004-02-01T00:00:00.000Z" [21,] "2004-03-01T00:00:00.000Z" [22,] "2004-04-01T00:00:00.000Z" [23,] "2004-05-01T00:00:00.000Z" [24,] "2004-06-01T00:00:00.000Z" [25,] "2004-07-01T00:00:00.000Z" [26,] "2004-08-01T00:00:00.000Z" [27,] "2004-09-01T00:00:00.000Z" [28,] "2004-10-01T00:00:00.000Z" [29,] "2004-11-01T00:00:00.000Z" [30,] "2004-12-01T00:00:00.000Z" [31,] "2005-01-01T00:00:00.000Z" [32,] "2005-02-01T00:00:00.000Z" [33,] "2005-03-01T00:00:00.000Z" [34,] "2005-04-01T00:00:00.000Z" [35,] "2005-05-01T00:00:00.000Z" [36,] "2005-06-01T00:00:00.000Z" [37,] "2005-07-01T00:00:00.000Z" [38,] "2005-08-01T00:00:00.000Z" [39,] "2005-09-01T00:00:00.000Z" [40,] "2005-10-01T00:00:00.000Z" [41,] "2005-11-01T00:00:00.000Z" [42,] "2005-12-01T00:00:00.000Z" [43,] "2006-01-01T00:00:00.000Z" [44,] "2006-02-01T00:00:00.000Z" [45,] "2006-03-01T00:00:00.000Z" [46,] "2006-04-01T00:00:00.000Z" [47,] "2006-05-01T00:00:00.000Z" [48,] "2006-06-01T00:00:00.000Z" [49,] "2006-07-01T00:00:00.000Z" [50,] "2006-08-01T00:00:00.000Z" [51,] "2006-09-01T00:00:00.000Z" [52,] "2006-10-01T00:00:00.000Z" [53,] "2006-11-01T00:00:00.000Z" [54,] "2006-12-01T00:00:00.000Z" [55,] "2007-01-01T00:00:00.000Z" [56,] "2007-02-01T00:00:00.000Z" [57,] "2007-03-01T00:00:00.000Z" [58,] "2007-04-01T00:00:00.000Z" [59,] "2007-05-01T00:00:00.000Z" [60,] "2007-06-01T00:00:00.000Z" [61,] "2007-07-01T00:00:00.000Z" [62,] "2007-08-01T00:00:00.000Z" [63,] "2007-09-01T00:00:00.000Z" [64,] "2007-10-01T00:00:00.000Z" [65,] "2007-11-01T00:00:00.000Z" [66,] "2007-12-01T00:00:00.000Z" [67,] "2008-01-01T00:00:00.000Z" [68,] "2008-02-01T00:00:00.000Z" [69,] "2008-03-01T00:00:00.000Z" [70,] "2008-04-01T00:00:00.000Z" [71,] "2008-05-01T00:00:00.000Z" [72,] "2008-06-01T00:00:00.000Z" [73,] "2008-07-01T00:00:00.000Z" [74,] "2008-08-01T00:00:00.000Z" [75,] "2008-09-01T00:00:00.000Z" [76,] "2008-10-01T00:00:00.000Z" [77,] "2008-11-01T00:00:00.000Z" [78,] "2008-12-01T00:00:00.000Z" [79,] "2009-01-01T00:00:00.000Z" [80,] "2009-02-01T00:00:00.000Z" [81,] "2009-03-01T00:00:00.000Z" [82,] "2009-04-01T00:00:00.000Z" [83,] "2009-05-01T00:00:00.000Z" [84,] "2009-06-01T00:00:00.000Z" [85,] "2009-07-01T00:00:00.000Z" [86,] "2009-08-01T00:00:00.000Z" [87,] "2009-09-01T00:00:00.000Z" [88,] "2009-10-01T00:00:00.000Z" [89,] "2009-11-01T00:00:00.000Z" [90,] "2009-12-01T00:00:00.000Z" [91,] "2010-01-01T00:00:00.000Z" [92,] "2010-02-01T00:00:00.000Z" [93,] "2010-03-01T00:00:00.000Z" [94,] "2010-04-01T00:00:00.000Z" [95,] "2010-05-01T00:00:00.000Z" [96,] "2010-06-01T00:00:00.000Z" [97,] "2010-07-01T00:00:00.000Z" [98,] "2010-08-01T00:00:00.000Z" [99,] "2010-09-01T00:00:00.000Z" [100,] "2010-10-01T00:00:00.000Z" [101,] "2010-11-01T00:00:00.000Z" [102,] "2010-12-01T00:00:00.000Z" [103,] "2011-01-01T00:00:00.000Z" [104,] "2011-02-01T00:00:00.000Z" [105,] "2011-03-01T00:00:00.000Z" [106,] "2011-04-01T00:00:00.000Z" [107,] "2011-05-01T00:00:00.000Z" [108,] "2011-06-01T00:00:00.000Z" [109,] "2011-07-01T00:00:00.000Z" [110,] "2011-08-01T00:00:00.000Z" [111,] "2011-09-01T00:00:00.000Z" [112,] "2011-10-01T00:00:00.000Z" [113,] "2011-11-01T00:00:00.000Z" [114,] "2011-12-01T00:00:00.000Z" [115,] "2012-01-01T00:00:00.000Z" [116,] "2012-02-01T00:00:00.000Z" [117,] "2012-03-01T00:00:00.000Z" [118,] "2012-04-01T00:00:00.000Z" [119,] "2012-05-01T00:00:00.000Z" [120,] "2012-06-01T00:00:00.000Z" [121,] "2012-07-01T00:00:00.000Z" [122,] "2012-08-01T00:00:00.000Z" [123,] "2012-09-01T00:00:00.000Z" [124,] "2012-10-01T00:00:00.000Z" [125,] "2012-11-01T00:00:00.000Z" [126,] "2012-12-01T00:00:00.000Z" [127,] "2013-01-01T00:00:00.000Z" [128,] "2013-02-01T00:00:00.000Z" [129,] "2013-03-01T00:00:00.000Z" [130,] "2013-04-01T00:00:00.000Z" [131,] "2013-05-01T00:00:00.000Z" [132,] "2013-06-01T00:00:00.000Z" [133,] "2013-07-01T00:00:00.000Z" [134,] "2013-08-01T00:00:00.000Z" [135,] "2013-09-01T00:00:00.000Z" [136,] "2013-10-01T00:00:00.000Z" [137,] "2013-11-01T00:00:00.000Z" [138,] "2013-12-01T00:00:00.000Z" [139,] "2014-01-01T00:00:00.000Z" [140,] "2014-02-01T00:00:00.000Z" [141,] "2014-03-01T00:00:00.000Z" [142,] "2014-04-01T00:00:00.000Z" [143,] "2014-05-01T00:00:00.000Z" [144,] "2014-06-01T00:00:00.000Z" [145,] "2014-07-01T00:00:00.000Z" [146,] "2014-08-01T00:00:00.000Z" [147,] "2014-09-01T00:00:00.000Z" [148,] "2014-10-01T00:00:00.000Z" [149,] "2014-11-01T00:00:00.000Z" [150,] "2014-12-01T00:00:00.000Z" [151,] "2015-01-01T00:00:00.000Z" [152,] "2015-02-01T00:00:00.000Z" [153,] "2015-03-01T00:00:00.000Z" [154,] "2015-04-01T00:00:00.000Z" [155,] "2015-05-01T00:00:00.000Z"
Get coverage's time axis's domain (as this coverage is 3D with time axis)¶
In [9]:
temporal_dim <- NULL ## Get coverage's time axis's domain (as this coverage is 3D with time axis)
for (dim in chla_dims) {
if (dim$type == "temporal") {
temporal_dim <- dim
}
}
if (!is.null(temporal_dim)) {
# Get the first and last time coefficients
time_coefficients <- temporal_dim$coefficients
first_time <- time_coefficients[1]
last_time <- time_coefficients[length(time_coefficients)]
print(paste("First time:", first_time))
print(paste("Last time:", last_time))
} else {
print("No temporal dimension found.")
}
[1] "First time: 2002-07-01T00:00:00.000Z" [1] "Last time: 2015-05-01T00:00:00.000Z"
Get a subset coverage¶
By slicing on time axis, trimming on Lat and Long axes, downloading as TIFF file, and displaying the result directly
In [28]:
data <- chla$getCoverage(
bbox = ows4R::OWSUtils$toBBOX(-122.1420, 122.2185 , -81.7242, 81.7825),
time = "2002-09-01T00:00:00.000Z",
filename = "myfile.tif"
) # Get a subset coverage by slicing on time axis, trimming on Lat and Long axes, and download as TIFF file
r <- raster::raster(file.path("~", "Downloads", "myfile.tiff"))
#Note: It was assumed that the file gets downloaded in the Downloads directory
raster::plot(r) # Display result directly
Tip: There is also an option to download multiple coverages in a stack, by slightly changing a previous function, the way as in the cell below:
In [24]:
cov_stack <- chla$getCoverageStack(
bbox = ows4R::OWSUtils$toBBOX(-122.1420,122.2185 , -81.7242, 81.7825),
time = tail(chla_dims[[1]]$coefficients, 5),
filename_handler = WCSCoverageFilenameHandler
) #Load stack of subset coverages
[ows4R][INFO] WCSCoverageSummary - Fetching coverage stack with 'temporal' dimension [ows4R][INFO] WCSCoverageSummary - 2015-01-01T00:00:00.000Z [ows4R][WARN] WCSCoverageSummary - Coverage without vertical dimension: 'elevation' argument is ignored <GMLEnvelope> ....|-- lowerCorner: "2002-07-01T00:00:00.000Z" -81.7242 -122.142 ....|-- upperCorner: "2015-05-01T00:00:00.000Z" 81.7825 122.2185[ows4R][INFO] WCSGetCoverage - Fetching https://ows.rasdaman.org/rasdaman/ows?service=WCS&version=2.1.0&coverageId=AverageChloroColorScaled&subset=ansi(%222015-01-01T00:00:00.000Z%22)&subset=Lat(-81.7242,81.7825)&subset=Lon(-122.142,122.2185)&format=image/tiff&request=GetCoverage Downloading: 760 kB [ows4R][INFO] WCSCoverageSummary - 2015-02-01T00:00:00.000Z [ows4R][WARN] WCSCoverageSummary - Coverage without vertical dimension: 'elevation' argument is ignored <GMLEnvelope> ....|-- lowerCorner: "2002-07-01T00:00:00.000Z" -81.7242 -122.142 ....|-- upperCorner: "2015-05-01T00:00:00.000Z" 81.7825 122.2185[ows4R][INFO] WCSGetCoverage - Fetching https://ows.rasdaman.org/rasdaman/ows?service=WCS&version=2.1.0&coverageId=AverageChloroColorScaled&subset=ansi(%222015-02-01T00:00:00.000Z%22)&subset=Lat(-81.7242,81.7825)&subset=Lon(-122.142,122.2185)&format=image/tiff&request=GetCoverage Downloading: 760 kB [ows4R][INFO] WCSCoverageSummary - 2015-03-01T00:00:00.000Z [ows4R][WARN] WCSCoverageSummary - Coverage without vertical dimension: 'elevation' argument is ignored <GMLEnvelope> ....|-- lowerCorner: "2002-07-01T00:00:00.000Z" -81.7242 -122.142 ....|-- upperCorner: "2015-05-01T00:00:00.000Z" 81.7825 122.2185[ows4R][INFO] WCSGetCoverage - Fetching https://ows.rasdaman.org/rasdaman/ows?service=WCS&version=2.1.0&coverageId=AverageChloroColorScaled&subset=ansi(%222015-03-01T00:00:00.000Z%22)&subset=Lat(-81.7242,81.7825)&subset=Lon(-122.142,122.2185)&format=image/tiff&request=GetCoverage Downloading: 760 kB [ows4R][INFO] WCSCoverageSummary - 2015-04-01T00:00:00.000Z [ows4R][WARN] WCSCoverageSummary - Coverage without vertical dimension: 'elevation' argument is ignored <GMLEnvelope> ....|-- lowerCorner: "2002-07-01T00:00:00.000Z" -81.7242 -122.142 ....|-- upperCorner: "2015-05-01T00:00:00.000Z" 81.7825 122.2185[ows4R][INFO] WCSGetCoverage - Fetching https://ows.rasdaman.org/rasdaman/ows?service=WCS&version=2.1.0&coverageId=AverageChloroColorScaled&subset=ansi(%222015-04-01T00:00:00.000Z%22)&subset=Lat(-81.7242,81.7825)&subset=Lon(-122.142,122.2185)&format=image/tiff&request=GetCoverage Downloading: 760 kB [ows4R][INFO] WCSCoverageSummary - 2015-05-01T00:00:00.000Z [ows4R][WARN] WCSCoverageSummary - Coverage without vertical dimension: 'elevation' argument is ignored <GMLEnvelope> ....|-- lowerCorner: "2002-07-01T00:00:00.000Z" -81.7242 -122.142 ....|-- upperCorner: "2015-05-01T00:00:00.000Z" 81.7825 122.2185[ows4R][INFO] WCSGetCoverage - Fetching https://ows.rasdaman.org/rasdaman/ows?service=WCS&version=2.1.0&coverageId=AverageChloroColorScaled&subset=ansi(%222015-05-01T00:00:00.000Z%22)&subset=Lat(-81.7242,81.7825)&subset=Lon(-122.142,122.2185)&format=image/tiff&request=GetCoverage Downloading: 760 kB
OGC Web Coverage Processing Cervices (WCPS) Requests¶
The OGC Web Coverage Processing Service (WCPS) is an extension of the OGC WCS standard that allows for the extraction, analysis, and processing of multi-dimensional coverages. WCPS is designed with syntax similar to the XQuery language, and it establishes a protocol for sending a query string to a server and receiving a set of coverages as the result of the server's processing.
The processing expression is applied to each coverage specified in the given list (coverageList) that satisfies the optional boolean expression. The query references each coverage using the corresponding identifier variableName in the processing expression. The processing expression can include sub-expressions that return scalars or encoded marrays, which operate on both the data and metadata of a coverage. More information can be found at http://tutorial.rasdaman.org/rasdaman-and-ogc-ws-tutorial/#ogc-web-services-web-coverage-processing-service.
To execute a query, use the following request (ProcessCoverages is an extension of the WCS request):
https://inspire.rasdaman.org/rasdaman/ows?service=WCS&version=2.0.1&request=ProcessCoverages&query=
The query parameter should contain a WCPS query with valid syntax for the server to process. An example of subsetting a 3D coverage to a 2D coverage and encoding the result in image/png is provided in the previous demo for the WCS GetCoverage request.
In [ ]:
# NOTE: for GET requests, WCPS query cannot contain new lines character (below query is just for brevity)
for $c in (AvgTemperatureColorScaled) return encode(
$c[Lat(-81.7242:81.7825),
Lon(-122.1420:122.2185),
ansi("2002-09-01T00:00:00.000Z")],
"png")
WCPS provides the capability for significantly more intricate processing requests than what can be achieved with WCS. For instance, it allows for time-series processing to determine the difference between two specific date-time slices ("2002-09-01T00:00:00.000Z" (reference) and "2009-05-01T00:00:00.000Z") of the coverage AvgTemperatureColorScaled.
**NOTE:** It should be noted that when using WCPS queries that contain special characters such as '[', ']', '{', '}', it is necessary to use POST requests.
Some of the libraries used are not available in CRAN R packages, have been found on GitHub, such as myCubeR (https://github.com/ARPASMR/myCubeR, and have been installed using devtools package's "install_github()" function.
Integrating a WCPS query and displaying the results¶
Slicing, trimming, and downloading a TIFF file¶
myCubeR::WPCS_query() function send a query to a server, and receives the results. The variables needed are the query itself, format, name of the file where we put the retreived data, and URL to which the query is sent. In fact you can send any type of an executable query.
In [22]:
query='for $c in (AvgTemperatureColorScaled) return encode(
($c[Lat(-81.7242:81.7825), Lon(-122.1420:122.2185), ansi("2002-09-01T00:00:00.000Z")]
- $c[Lat(-81.7242:81.7825), Lon(-122.1420:122.2185), ansi("2009-05-01T00:00:00.000Z")])
* 200, "tiff")' # Query for a TIFF format
raster_tot=myCubeR::WPCS_query(proper_query=query, FORMAT="image/tiff", filename="AvgTemperatureColorScaled.tiff", query_url="https://inspire.rasdaman.org/rasdaman/ows")
# Send a WCPS query with special characters to server and download the result
file_path <- file.path("~", "Downloads", "AvgTemperatureColorScaled.tiff")
ra <- raster::raster(file_path)
#Note: It was assumed that the file gets downloaded in the Downloads directory
raster::plot(ra)
# Display a TIFF file directly
Tip: The data can also be downloaded and displayed in a PNG or JPEG formats. Here we display the same query with a PNG encoding
In [18]:
query='for $c in (AvgTemperatureColorScaled) return encode(
($c[Lat(-81.7242:81.7825), Lon(-122.1420:122.2185), ansi("2002-09-01T00:00:00.000Z")]
- $c[Lat(-81.7242:81.7825), Lon(-122.1420:122.2185), ansi("2009-05-01T00:00:00.000Z")])
* 200, "png")'
raster_top=myCubeR::WPCS_query(proper_query=query, FORMAT="image/png", filename="AvgTemperatureColorScaled.png", query_url="https://inspire.rasdaman.org/rasdaman/ows")
# Send a WCPS query with special characters to server and download the result
file_path = file.path("~", "Downloads", "AvgTemperatureColorScaled.png")
img <- png::readPNG(file_path)
#Note: It was assumed that the file gets downloaded in the Downloads directory
grid::grid.raster(img)