Multidimensional Visualization and Processing of Big Open Urban Geospatial Data on the Web

: The focus of this research is addressing a subset of the geovisualization (i.e., geographic visualization) challenges identiﬁed in the literature, namely multidimensional vector and raster geospatial data visualization. Moreover, the work implements an approach for multidimensional raster geospatial data processing. The results of this research are provided through a geoportal comprised of multiple applications that are related to 3D visualization of cities, ground deformation, land use and land cover and mobility. In a subset of the applications, the datasets handled are considered to be large in volume. The geospatial data were visualized on dynamic and interactive virtual globes to enable visual exploration. The geoportal is available on the web to enable cross-platform access to it. Furthermore, the geoportal was developed employing open standards, free and open source software (FOSS) and open data, most importantly to ensure interoperability and reduce the barriers to access it. The geoportal brings together various datasets, different both in terms of context and format employing numerous technologies. As a result, the existing web technologies for geovisualization and geospatial data processing were examined and exemplary and innovative software was developed to extend the state of the art.


Introduction
Geographic visualization, often shortened as geovisualization, was born as a result of the tools emerged in the digital era, transforming cartography. Geovisualization integrates approaches from visualization in scientific computing, cartography, image analysis, information visualization, exploratory data analysis (EDA) and geographic information systems (GIS) to provide theory, methods and tools for visual exploration, analysis, synthesis and presentation of geospatial data [1]. Geovisualization techniques have extended the static 2D products and 3D models by introducing dynamic and interactive 3D and 4D data representation. Such representations are capable of incorporating augmented and virtual realities and geodatabases. Moreover, they are flexible in application, platform, scale and content and the data that they represent are possibly multivariate [2]. The most distinct feature of geovisualization compared to traditional cartography is its interactivity, which enables exploration and thus making inferences [3].
The potential of using virtual globes for visualization (i.e., 3D and realistic geospatial data visualization) is a paradigm that first emerged with the development of NASA WorldWind in 2003. Currently, there are inadequate applications designed to utilize the most recent free and open source technologies that respond to this paradigm for the realization of the Digital Earth [4][5][6][7][8][9]. Furthermore, not restricted to the themes of the applications in this work, it is envisaged that a single system that brings together various important aspects of urban areas will be desirable and the datasets employed will likely to have multidimensional nature, as there is a trend towards both multidimensional vector data visualization (i.e., creating digital twins of cities) and having the time dimension in both visualization and processing of raster data as a result of the continuous production of EO data. The authors could have chosen to have all the applications in a single virtual globe application, but they created sections according to the themes to facilitate the navigation. Having a single entry point to a vast amount of data under various themes through the geoportal aims to both communicate the outputs of the research carried out in the URBAN GEO BIG DATA project and by doing so reduce the gap between science and society and aid scientists, public and private sector administrators and citizens in decision-making regarding various real-world phenomena through dynamic and interactive visual interfaces that provide an effective way for human reasoning. Overall, this work aims to contribute to the vision of Digital Earth.

An Overview of Technologies
The geospatial visual interface of the geoportal of the URBAN GEO BIG DATA project is composed of multiple dynamic and interactive virtual globes for different applications. The geoportal, besides visualization, enables the query and processing of certain datasets to incite human reasoning and aid in effective decision-making further. As research groups in different fields produced the datasets, bringing them together in the same geoportal poses technical challenges that were solved to produce the geoportal. More specifically, the datasets are multidimensional, challenging to handle due to their large sizes and both in vector and raster format.
The geoportal is available on the web so that access to it is independent of users' operating system. Moreover, open standards were adopted to ensure the interoperability of the geoportal. Furthermore, free and open source software (FOSS) and open data were used for developing the geoportal most importantly to reduce the barriers to access it.
The open standards and FOSS that were used for tackling the challenges presented in the previous section through developing the geoportal of the URBAN GEO BIG DATA project are described in this section. The standards and software cover a wide range of functionalities, primarily storage, visualization, query and processing of big multidimensional geospatial data in vector or raster format. Their strategies for dealing with such data that were utilized include tiling and indexing to access data faster and parallel processing. The standards and software were chosen in a way to ensure interoperability of and cross-platform access to the geoportal. The overview of the project can be obtained from Figure 1. Software architecture of the URBAN GEOmatics for Bulk Information Generation, Data Assessment and Technology Awareness (URBAN GEO BIG DATA) project that is composed of a geoportal for geospatial data visualization, query and processing, a discovery service for enabling the search of geospatial data based on various parameters and data servers that host the geospatial data used by the geoportal and the discovery service.

Open Standards
The application programming interfaces (APIs) employed for geovisualization that are described in the Section 2.2 use Web Graphics Library (WebGL), a cross-platform and royalty-free web standard for a low-level graphics API, designed and maintained by the non-profit Khronos Group. The implementation of the standard is a JavaScript API for rendering interactive 2D and 3D graphics within a web browser without a plug-in as the API closely conforms to Open Graphics Library for Embedded Systems (OpenGL ES) that can be used in Hypertext Markup Language 5 (HTML5) canvas elements.
The geospatial data are retrieved by the client from several servers using the open web standards of Open Geospatial Consortium (OGC). The OGC standards deliver or process the raw geospatial data or deliver the portrayal of the raw geospatial data stored on a server or multiple servers. The OGC standards used are Web Map Service (WMS), Web Map Tile Service (WMTS), Web Feature Service (WFS), Web Coverage Service (WCS) and WCPS.
WMS produces maps of georeferenced data dynamically from geospatial data. A map is a portrayal of geospatial data as a digital image suitable for display on a computer screen, as a result, it is not the data. Styled Layer Descriptor (SLD) defines how a WMS can be extended to allow user-defined symbolization and colouring of geographic features and coverage data using Extensible Markup Language (XML) encoding. SLD is an OGC specification. WMTS provides a complementary approach to WMS. WMS generates a map dynamically, possibly using an SLD document for rendering with a custom style. WMTS trades the flexibility of custom map rendering for the scalability possible by serving predefined image tiles. Serving static data allows the implementation of a WMTS using a web server that simply returns existing files and the use of standard network mechanisms for scalability such as distributed cache systems. On the contrary, WFS offers access to geographic data at the feature and feature property level. It allows clients to retrieve, modify, replace and delete the data they are seeking, rather than a file that contains the data they are seeking.
The WCS has a mandatory Core, around which extensions can be implemented. WCS Core offers coverage access, subsetting and format encoding. Subsetting a coverage allows trimming and slicing. Trimming returns a coverage with the same number of dimensions as the input coverage. On the other hand, slicing returns a coverage with a reduced number of dimensions compared to the input coverage (see Figure 2). The list of WCS Core extensions can be found in [33]. The one that is relevant to this is work is WCPS. WCPS defines a protocol-independent language for the extraction, processing and analysis of multidimensional coverages representing sensor, image or statistics data. It allows coverage analytics through appending a single parameter query to the WCS request, which contains a WCPS query string. Concisely, WCPS allows: • Subsetting (downloading a subset of a coverage by trimming or slicing); • Range subsetting (extracting a band or bands of a coverage); • Condensing (consolidating cell values of a coverage along selected axes to a scalar value based on a condensing operation, such as calculating minimum, maximum, average or sum of the cell values); • Constructing a coverage (creating a new coverage on the fly and filling it with values resulting from a processing expression evaluation); • Applying induced operations (using a unary or binary function, that may include arithmetic, comparison, Boolean, trigonometric and logarithmic operations and case distinction that works on a single cell and applying it to all the cells of a coverage simultaneously).
As coverages are a superset of datacubes, the standards related to coverages are used to implement services for accessing and processing datacubes. Ref. [33] defined a datacube as: A datacube is a massive multidimensional array, also called 'raster data' or 'gridded data'; 'massive' entails [...] sizes significantly beyond the main memory resources of the server hardware-otherwise, processing can be done satisfactorily with existing array tools like MATLAB or R.
The datacube concept heralds increased access to EO data, including additionally the time axis. Datacubes simplify access to multi-temporal geospatial data by constituting a 3D array (x/y/t) for raster time series. Moreover, datacubes may constitute 1D sensor time series, 2D imagery, 3D (x/y/z) subsurface voxel data, 4D (x/y/z/t) climate and ocean datacubes and even 5D atmospheric data with two time dimensions [33]. In this work, 2D imagery and 3D array (x/y/t) for raster time series were used. As a result, hereafter, whenever the term coverage is used in relation to the WCPS standard, the data are in raster format.
Moreover, in terms of open data standards, City Geography Markup Language (CityGML), JavaScript Object Notation (JSON), GeoJSON and Geographic Tagged Image File Format (GeoTIFF) were used. CityGML is an OGC encoding standard and represents an open data model and XML-based format for virtual 3D city models. Furthermore, GeoJSON is an open standard geospatial data format based on JSON.

Free and Open Source Software
There are various APIs for creating virtual globes on the web. Among them, the two most prominent ones are NASA Web WorldWind and CesiumJS. Both APIs offer high-performance, dynamic, interactive and cross-platform visualization and exploration of a virtual globe with terrain and 2D and 3D geospatial data layers on it. Both utilize WebGL, JavaScript, HTML5 and Cascading Style Sheets Level 3 (CSS3).
Today, many proprietary and open source tools support CityGML. An exhaustive list regarding these tools can be found in the wiki of CityGML (http://www.citygmlwiki.org/). The 3D City Database (3DCityDB) is one of these tools. It is a geodatabase for storing and managing virtual 3D city models on top of a standard spatial relational database (Oracle Spatial or Locator or PostgreSQL with PostGIS). The database schema fully implements the CityGML 1.0.0 and 2.0.0 standards. The database model contains semantically rich, hierarchically structured, multi-scale city objects facilitating complex GIS modelling and analysis tasks, beyond visualization. 3DCityDB Importer/Exporter is a Java-based application for importing data in CityGML format and exporting tiled data in Keyhole Markup Language (KML), Collaborative Design Activity (COLLADA), GL Transmission Format (glTF) or CityGML formats. 3DCityDB-Web-Map-Client extends CesiumJS for 3D visualization and interactive exploration of arbitrarily large semantic virtual 3D city models on the web. It allows visualizing large and tiled data in KML, COLLADA, glTF formats obtained using the 3DCityDB Importer/Exporter along with terrain and vector and raster layers.
GeoServer is a Java-based FOSS server for sharing geospatial data using open standards. It implements WMS, Transactional WFS and WCS. Additionally, WPS, WMTS and Catalogue Service for the web (CSW) implementations are provided as extensions. The clients that request the geospatial data published on GeoServer are typically web browsers and desktop GIS software. GeoServer has a management interface available on a web browser that connects to the data sources at the back-end. Moreover, GeoWebCache (GWC) is integrated into GeoServer. GWC is a Java-based free and open source web application for caching map tiles coming from various sources. rasdaman (raster data manager) has been developed for more than two decades into a cross-domain datacube engine over a series of projects [34][35][36][37][38][39][40]. The raster query language of rasdaman, rasql extends Structured Query Language (SQL) with declarative nD array operators [41]. A separate layer adds geo semantics, such as information about regular and irregular grids and coordinate reference systems (CRSs), by implementing WMS 1.3, WCS 2.0 and WCPS 1.0. For OGC and INSPIRE WCS and OGC WCPS rasdaman is the reference implementation. rasdaman is the only available implementation for OGC WCPS.

Results
In the following sections, each application is explained in detail with the datasets they employed, their software development details and functionalities. When the background or literature diverges from the one given in the first two sections, they are given in a concise manner.

OSM Data Visualization
OpenStreetMap (OSM) was founded in 2004 and has become a collaborative project that created a free and editable map of the world. OSM is one of the most prominent VGI projects. VGI, coined by [42], is a specific case of the more general phenomenon, user-generated content (UGC) from Web 2.0. VGI is geospatial UGC, i.e., geospatial data collected by individuals voluntarily [43]. OSM database contains crowdsourced geospatial data that are open, licensed under Open Database License (ODbL) [44].
Since 2008, people more and more often have collected data that can be used for creating 3D objects, such as height and roof geometry data [45]. Being able to constitute 3D objects from VGI sources is important, not only because we live in a 3D world, but also because 3D data allow the development and provision of many applications. The applications include providing a virtual 3D model of a city for demonstrating future city development plans to the public [46], visibility analysis [47], decision support for emergency [48][49][50], visualizing topological relations [51] and determining the environmental quality of public spaces [52]. Various applications use OSM data for building 3D objects. A comprehensive list is available in the wiki of OSM (https://wiki.openstreetmap.org/wiki/3D). Moreover, extensive research has been carried out related to this topic [45,[53][54][55][56].
The 3D OSM Plugin API was developed for NASA Web WorldWind in the Google Summer of Code (GSoC) program in 2017. The source code of the API is available on GitHub (https://github. com/kilsedar/3dosm) with MIT License. The API provides a way to visualize OSM data in 2D and 3D on a virtual globe created using NASA Web WorldWind. Only buildings are visualized in 3D; the rest of the features are visualized in 2D. The spatial quality of OSM buildings was checked by comparing the OSM data of the Lombardy region of Italy whose capital is Milan with the region's authoritative dataset. As a result, it was found that the positional accuracy of the OSM buildings is comparable with the quality of the authoritative dataset at scale 1:5000, as a result, they are suitable to use for the visualization developed in this work [57].
The 3D OSM Plugin API fetches the OSM data based on a tag and the bounding box of the area for which the user wants to retrieve elements. Overpass API was used for fetching the data. Overpass API is a read-only API that serves a selected part of the OSM database. A client sends a query to the API and gets the data that correspond to the query. It has two query languages: Overpass XML and Overpass Query Language (Overpass QL). The Overpass QL was used to construct the queries. Overpass API returns the data in JSON format. The data in JSON format were converted to GeoJSON format using the osmtogeojson API (https://github.com/tyrasd/osmtogeojson). Besides getting data using the Overpass API, it is possible to use a file in GeoJSON format or write the data in GeoJSON format inside JavaScript code. Only the GeoJSON data obtained in the way that the 3D OSM Plugin API does were tested and the results were as expected.
The height of the buildings that are visualized in 3D can be set using the OSM database, more specifically, the tags of features. First, the height key (https://wiki.openstreetmap.org/wiki/Key: height) that describes the height of a feature is used to set the height of a building, if a value is assigned to it for an element with the building key (https://wiki.openstreetmap.org/wiki/Key:building). If a value for such a height key is not defined, the building:levels key (https://wiki.openstreetmap.org/ wiki/Key:building:levels) that describes the number of above-ground levels of a building or part of a building is used to set the height of a building, if a value is assigned to it. Each level is assumed to be 3 m. If a value is not assigned to either key, the building is assumed to have five levels, i.e., be 15 m long. A comprehensive list of OSM keys that can be used to render 3D buildings is available in the wiki of OSM (https://wiki.openstreetmap.org/wiki/Simple_3D_buildings). As these tags are not present most of the time, the API lets the user use a file in GeoJSON format obtained in the way described above, with an additional attribute that denotes the building heights. According to the publication [45], less than 1.5% of the elements with building key in the OSM database had the height key in November 2011. Moreover, it is possible to set an arbitrary value for all the buildings. It is also possible to visualize the buildings in 2D.
The style of the OSM elements can be set using the API that employs two classes of NASA Web WorldWind: PlacemarkAttributes and ShapeAttributes. PlacemarkAttributes class is used for Point and MultiPoint Geometry types of GeoJSON. ShapeAttributes class is used for the rest of the Geometry types of GeoJSON and triangle meshes. Using the PlacemarkAttributes class, it is possible to set the image and label of placemarks, among others. Using the ShapeAttributes class, it is possible to set the interior and outline colour, outline width, whether to draw an outline, among others. The user can choose to set varying colours for the 3D visualization of buildings. The colours depend on the height of the buildings in the dataset and the colour and thresholds defined by the user. As the height of the buildings increases, the red component in the colour in RGB gets a higher value. Moreover, the thresholds that are ordered hold the values for building heights, and single colour is assigned to all the buildings with a height between two consecutive thresholds. In this way, a 3D heatmap that represents building heights is generated.
The authors participated in NASA Europa Challenge 2017 with an application that uses the 3D OSM Plugin API. Figure 3 displays screenshots from the application and exemplifies the capabilities of the API. The elements that are in the bounding box drawn on the virtual globe by the user can be retrieved from the OSM database. Moreover, the tag and colour of the elements can be defined by the user using the graphical user interface (GUI) of the application. Figure 3a shows the amenities in Istanbul displayed using placemarks. The nodes were retrieved from the OSM database using the amenity=yes tag (https://wiki.openstreetmap.org/wiki/Key:amenity). Figure 3b shows the pathways that are used exclusively or mainly by pedestrians in Milan displayed using lines and polygons. The ways were retrieved from the OSM database using the highway=footway tag (https: //wiki.openstreetmap.org/wiki/Tag:highway=footway). Figure 3c shows the forests and woodlands in Helsinki displayed using polygons and multipolygons. The ways and relations were retrieved from the OSM database using the landuse=forest tag (https://wiki.openstreetmap.org/wiki/Tag: landuse=forest). Figure 3d shows the buildings in New York displayed using triangle meshes based on polygons and multipolygons. The ways and relations were retrieved from the OSM database using the building=yes tag. This visualization uses the values of the keys height and building: levels associated with the elements that represent a building to set the height of the triangle meshes. This visualization clearly shows the varying height of the buildings, which is because for most of the elements that represent a building in New York the height or building:levels keys with an assigned value exist.
(a) Amenities displayed using placemarks  However, as most of the elements that have the building key do not have the height or building:levels key, when possible, an alternative way was used to set the building heights. First, LIDAR data that represent digital terrain model (DTM) and digital surface model (DSM) were used to extract the height of the buildings in Milan. LIDAR data were received from the Ministry for Environment, Land and Sea Protection of Italy using the contact information available at the national geoportal (http://www.pcn.minambiente.it/mattm/en/data-distribution-service-pst/). Then, they were converted to virtual datasets (VRTs) using the gdalbuildvrt program of Geospatial Data Abstraction Library (GDAL). Then, the VRTs of DTM and DSM were imported into GRASS GIS. The procedure involves subtracting DTM from DSM to get the height of the objects on the terrain using the r.mapcalc module of GRASS GIS. The OSM data that contain the building footprints in Milan were downloaded in GeoJSON format using Mapzen Metro Extracts and imported into GRASS GIS. Using the v.rast.stats module of GRASS GIS, some statistics were calculated for the footprints of the buildings. The calculated statistics were the minimum, maximum, average and median of the values of pixels inside each footprint. Among these statistics, the median was used to discard the outliers that may stem from the misalignment between the footprints in the OSM data and LIDAR measurements. Second, Urban Atlas Building Height 2012 dataset in GeoTIFF format provided by the Copernicus Programme (https://land.copernicus.eu/local/urban-atlas/building-height-2012) was used to extract the height of the buildings in Rome. The only difference between this procedure and the previous one is the removal of subtraction between DTM and DSM, as the raster data already represent the building heights. For the other three cities the URBAN GEO BIG DATA project focuses on (i.e., Padua, Turin and Naples), the height and building:levels keys associated with the elements that represent a building were used to set the building heights.

CityGML Data Visualization
Besides VGI (i.e., OSM data in GeoJSON format), 3D vector data visualization was achieved using the OGC standard for representing 3D vector data that pertain to urban areas, CityGML. The visualization of the data in CityGML format would ideally be achieved using the 3D Tiles standard, as it is designed for streaming and rendering massive 3D geospatial data and using an OGC standard is consistent with the goal of achieving interoperability. Moreover, CesiumJS supports the 3D Tiles format. It is important to note that there are other OGC standards for 3D geospatial data delivery, such as 3D Portrayal Service (3DPS) and Indexed 3D Scene Layers (I3S), yet they are not supported by CesiumJS. However, there is not FOSS for converting data in CityGML format to 3D Tiles format. As a result, 3DCityDB was used for the visualization of the data in CityGML format in the URBAN GEO BIG DATA project.
Virtual 3D city models in CityGML format have been created for many cities in the world, yet they are not available for the five cities that the URBAN GEO BIG DATA project focuses on. For this reason, the team from the University of Padua developed software called shp2city that converts data in Esri shapefile format to CityGML format [58]. In the URBAN GEO BIG DATA project, CityGML datasets that represent the buildings of Milan, Padua, Turin and Naples were generated using the shp2city software. Using 3DCityDB Importer/Exporter, the datasets were imported into a PostgreSQL database extended with PostGIS. Then, using the same software they were exported in KML, COLLADA, glTF formats tiled. The exported datasets were visualized using the 3DCityDB-Web-Map-Client. Since the buildings in the CityGML datasets have altitude values, it is possible to place them on the terrain. On the geoportal, it is also possible to simulate the sun, which enables to visualize shadows of the terrain and buildings at different times of the day and year. Moreover, six base maps were provided in the geoportal when CesiumJS was used to create the virtual globe, which are Bing Maps Aerial, Mapbox Satellite Streets, OSM, CARTO Dark, Stamen Terrain and Stamen Watercolor. The geoportal enables to switch between these six base maps.
Virtual 3D city models in CityGML format have been used for various applications, such as noise propagation simulation and mapping, energy-related assessments of buildings, indoor navigation, disaster management and homeland security [59,60]. Disaster management applications include flood simulations for assessing the flood risk and potential damage at a micro-scale [61][62][63][64][65]. In this work, virtual 3D city models in CityGML format were used for flood simulation. A semi-transparent polygon was placed on the ellipsoid surface of the virtual globe that users can extrude in meters using a slider in the GUI of the geoportal. It is possible to use the VRTheWorldTerrainProvider class of CesiumJS to generate the terrain geometry by tessellating the digital elevation model (DEM) that includes both land topography and bathymetry with a 90-m resolution for the entire globe provided by the VR-TheWorld Server (https://www.mak.com/products/terrain/vr-theworld-server). Instead of using the DEM provided by the VR-TheWorld Server, a 5-m DTM of Milan in GeoTIFF format was used to construct the terrain of the virtual globe to increase the accuracy of the flood simulation. This dataset is published by the Lombardy region (http://www.geoportale.regione.lombardia.it/en/home) as open data. The first step for using the DTM of Milan to construct the terrain of the virtual globe was creating terrain tiles in quantized-mesh-1.0 format (https://github.com/AnalyticalGraphicsInc/quantized-mesh). The tiles were created using the Cesium Terrain Builder (https://github.com/ahuarte47/cesiumterrain-builder). Then, the Cesium Terrain Server (https://github.com/geo-data/cesium-terrainserver) was used to host the tiles. The geoportal refers to this server to retrieve the tiles. Following this approach, it is not possible to simulate floods anywhere else in the world as the DTM is only for Milan unless additional local DEM or DTM is used to construct the terrain of the virtual globe. Finally, the flood risk map of Milan published by the Lombardy region as open data was stored on GeoServer and tiled using GWC. The ImageryLayer and WebMapTileServiceImageryProvider classes of CesiumJS were used by the client to make requests to the WMTS implementation of GeoServer and place the retrieved tiled imagery on the virtual globe. More information was given by [66]. Figure 4 summarizes the software used and their interactions to clarify the content given above. Figure 5 displays a flood simulation in Milan with the flood risk map and the visualization of the CityGML data that represent buildings.

Ground Deformation Visualization and Query
Differential Synthetic Aperture Radar Interferometry (DInSAR) is a well-established technique for the mapping and continuous monitoring of the areas on Earth that are subject to ground deformation [67]. National Research Council of Italy (CNR) Institute for Electromagnetic Sensing of the Environment (IREA) used the Small BAseline Subset (SBAS) technique, which is a DInSAR algorithm [68], to generate the mean deformation velocity maps and deformation time series for Milan, Padua, Turin, Naples and Rome. The deformation data give information on the changes of the Earth's surface with a direction and magnitude in a metric unit, which enables monitoring natural hazards and environmental alterations due to subsidence and agriculture, among others. This information facilitates determining the parts of the monitored areas that are susceptible to damage in case of natural hazards and the structures that are susceptible to damage or collapse as a result of or independent from natural hazards.
The SBAS algorithm was applied to the sequences of archived synthetic aperture radar (SAR) images collected by the European Remote Sensing (ERS) and Environmental Satellite (Envisat) satellites of ESA from 1992 to 2010. A maximum perpendicular baseline of 400 m and a maximum period of two years were used for the selection of the small baseline interferometric data pairs for the areas of interest in the five cities. Deformation data were produced at medium spatial resolution and for some parts of the residential and central regions of the areas of interest also at full spatial resolution.
The mean deformation velocity maps and deformation time series of the five cities are available on a Geoinformation Enabling ToolkIT starterkit R (GET-IT) (http://www.get-it.it/) installation (https: //ugbd.get-it.it/). GET-IT is developed at CNR IREA within the flagship project RITMARE (http: //www.ritmare.it/en/) starting from a widely known geospatial content management system (CMS) GeoNode [69]. GET-IT can be used to visualize and download geospatial data. Moreover, it can be used to read and download the metadata of the geospatial data. The visualization of the mean deformation velocity maps for the five cities is given in Figure 6. A mean deformation velocity map was retrieved by making a request to the GET-IT instance of the URBAN GEO BIG DATA project, using the WMS implementation of GeoServer. The ImageryLayer and WebMapServiceImageryProvider classes of CesiumJS were used by the client to make the request and place the retrieved image on the virtual globe. Each deformation point can be clicked on to display its cumulative deformation time series using a plot created employing a JavaScript API Plotly (https://github.com/plotly/plotly.js). The data used to create the plot were retrieved in JSON format by making a request to the GET-IT instance of the URBAN GEO BIG DATA project, using the WFS implementation of GeoServer. The style of the maps was defined using SLD encoding in GeoServer. The colour scheme for mean deformation velocity used in the URBAN GEO BIG DATA project is given in Figure 6. The visualization of the mean deformation velocity map of Rome and the plot of cumulative deformation time series of a deformation point is displayed in Figure 7. Plotting cumulative deformation time series is available for all the deformation points of the five cities. Furthermore, raster files that represent cumulative deformation were generated for almost every year and the five cities, so that an overview of deformation for eighteen years can be achieved by temporally animating raster images. The raster files in GeoTIFF format were derived from the deformation time series. Ref. [70] explained the method for creating the raster files, which have reduced spatial and temporal resolution; as a result, the deformation time series were summarized. The ImageMosaic data store of GeoServer was used for creating raster time series. An ImageMosaic was created for each of the five cities. The colour scheme for cumulative deformation used in the URBAN GEO BIG DATA project and the visualization of the raster deformation time series of Turin are given in Figure 8. The figures show the decreasing, increasing and periodic changes of deformation. An animation is available for each ImageMosaic. An ImageMosaic was temporally animated employing the animation and timeline widgets of CesiumJS created using the library's Animation and Timeline classes. The raster time series in ImageMosaic format comprises of multiple raster files, each having a timestamp. Thus, an animation is available for each ImageMosaic created in this work where each frame of the animation corresponds to a raster image with a timestamp. The animation can be initiated, paused and played forward and reverse using the animation widget. Moreover, using the timeline widget, it is possible to scroll through time manually. Animation and timeline widgets and the visualization of the raster file of 2005 from the raster deformation time series of Turin can be seen in Figure 9. The geoportal allows zooming to each of the five cities using a drop-down list. Once a city is selected using the list, either the mean deformation velocity map or the raster deformation time series of the selected city can be added to the virtual globe. As the deformation data have both negative and positive values, the colour schemes adhere to the diverging schemes where equal emphasis is put on the mid-range and extremes at both ends of the data range. In diverging schemes, the break in the middle is emphasized with light colours, and the low and high extremes are emphasized with dark colours that have contrasting hues [71]. Diverging schemes with eleven classes were chosen using ColorBrewer 2.0 (http://colorbrewer2.org) [72]. The mid-range value was replaced with 0 in the colour schemes to define different colours for negative and positive values. Moreover, the intervals for colours in the colour scheme for mean deformation velocity do not exceed the standard deviation of about 1 mm/year for mean deformation velocity. Likewise, the intervals for colours in the colour scheme for cumulative deformation do not exceed the standard deviation of about 5 mm for deformation [73]. Considering the standard deviation and the diverging nature of the data, the colour schemes given in Figures 6 and 8 were used. Both colour schemes were used for all the five cities to enable visual comparison among them. As the values for both mean deformation velocity and cumulative deformation concentrate around zero, the intervals become wider towards the low and high extremes.

LULC Visualization, Query and Processing
The growth of cities and as a result the increase of soil consumption [74] results in the decrease of life quality in the cities. One way it does so is by increasing the air temperature [75] which exacerbates health problems, especially for elders, during summertime. As a result, this work also focused on land use and land cover (LULC). EO-derived LULC datasets were visualized on a virtual globe on the web. Moreover, query and processing of the datasets can be initiated using the geoportal. The GlobeLand30 and GHS-BUILT that are available globally were cropped to the boundary of Italy. The size of the land consumption maps from ISPRA is around 224 GB, the size of the land cover map from ISPRA is around 28 GB, the size of GlobeLand30 is around 2.94 GB and the size of GHS-BUILT is around 6.23 GB. The land consumption maps from ISPRA are considered to be land use maps, and the rest are considered to be land cover maps. The land consumption maps from ISPRA and GHS-BUILT have two classes that indicate if the area is consumed or not and is built-up or not, respectively. The land cover map from ISPRA and GlobeLand30 have more than two classes.
The style of all the LULC datasets was defined using SLD encoding in GeoServer. The geoportal allows zooming to each of the five cities using a drop-down list. The land cover map from ISPRA of 2012 was stored using the GeoTIFF data store, instead of the ImageMosaic data store of GeoServer as it is available for a single year. After publishing the dataset on GeoServer, it was tiled using GWC. The ImageryLayer and WebMapTileServiceImageryProvider classes of CesiumJS were used by the client to make requests to the WMTS implementation of GeoServer and place the retrieved tiled imagery on the virtual globe.
As already mentioned in the previous section, the ImageMosaic data store of GeoServer was used for creating raster time series. An ImageMosaic was created for the land consumption maps from ISPRA, GlobeLand30 and GHS-BUILT. An animation is available for each ImageMosaic, which was created using the method described in the previous section. GlobeLand30 and GHS-BUILT that are available globally were cropped to the boundary of Italy as the URBAN GEO BIG DATA project focuses on Italy and the limited storage capacity. Figure 10 shows that animation allows detecting the changes of land cover over time visually as it is visible that built-up areas increased steadily in Rome from 1975 to 2014 according to GHS-BUILT. Furthermore, all the LULC datasets were imported into rasdaman so that they can be processed using WCPS. The datasets were imported using the web Coverage Service Transaction (WCS-T) implementation of rasdaman. As a result of the import of the sets of the datasets, a datacube was created that corresponds to each set of the datasets. Consequently, four datacubes were created. The datacube of the land cover map from ISPRA is 2D (x/y), while the rest of the datacubes are 3D (x/y/t).
Processing can be executed for both a set of coordinates and four sets of coordinates, i.e., a rectangle. When a user clicks on a pixel, the pixel coordinates of the canvas are translated into a set of coordinates in EPSG:3857 using CesiumJS. The processing, which represents a slicing operation, uses the set of coordinates and executes differently depending on the number of dimensions of the datacube. If the datacube is 2D, which is true for the land cover map from ISPRA, the processing returns the value of the pixel that is the land cover class. On the other hand, if the datacube is 3D, which is true for the land consumption maps from ISPRA, GlobeLand30 and GHS-BUILT, the processing returns the value of the pixel for all the times that the datasets are available, which corresponds to the land use or land cover classes. As a result, it is possible to get the change of land use or land cover classes over time for a pixel (see Figure 11). As mentioned earlier, the processing can also be executed for four sets of coordinates, besides a set of coordinates. The geoportal provides an interface for calculating the amount of change of a selected land use or land cover class inside a rectangle drawn by the user between two selected years. This operation is available only for 3D datacubes. The processing represents a combination of slicing and condensing operations. The four sets of pixel coordinates of the canvas that define the vertices of the rectangle are translated into four sets of coordinates in EPSG:3857 using CesiumJS.
The URLs executed (i.e., the requests made to the WCPS implementation of rasdaman) are http: //urbangeobigdata.como.polimi.it:8081/rasdaman/ows?query=for$cin(ispra_lc_2012_2015_2016_2017) returnencode(count($c[X(xMin:xMax),Y(yMin:yMax),ansi(firstYear)]=classificationCode),''text/csv'') and http://urbangeobigdata.como.polimi.it:8081/rasdaman/ows?query=for$cin(ispra_lc_2012_2015_ 2016_2017)returnencode(count($c[X(xMin:xMax),Y(yMin:yMax),ansi(secondYear)]=classificationCode) ,''text/csv'') for processing the datacube of the land consumption maps from ISPRA for four sets of coordinates. The xMin variable stands for minimum longitude, xMax variable stands for maximum longitude, yMin variable stands for minimum latitude and yMax variable stands for maximum latitude. These variables have the values of the four sets of coordinates in EPSG:3857 that correspond to the vertices of the rectangle. The value of the firstYear variable is the first year and the value of the secondYear variable is the second year the user chooses using the drop-down lists in the GUI. The classificationCode variable has the value of the pixel in the raster file that corresponds to the class the user chooses using the drop-down list in the GUI. The query parameter's value is a rasql query. For the other two datacubes, only the coverage id (i.e., ispra_lc_2012_2015_2016_2017) is different in the URLs. The requests calculate the number of pixels inside the rectangle drawn by the user that has the value of the classificationCode for the values of both the firstYear and secondYear. After getting the counts for both years, the percentage change is calculated and given to the user in a window created using the InfoBox class of CesiumJS (see Figure 12). VGI was visualized on the LULC raster imagery on the virtual globe. There are plenty of other applications for the visualization of VGI on a 2D or 3D map [76][77][78]. In this work, the VGI is collected by an application called Land Cover Collector that enables the collection of land cover data according to GlobeLand30 nomenclature.
The VGI is retrieved from a CouchDB database in JSON format. Each object in the JSON data that corresponds to a single submission of VGI is represented by an object created using the Entity class of CesiumJS. The position of an Entity is determined by the values of the longitude and latitude properties of the JSON object that it corresponds to. For each Entity, BillboardGraphics class of CesiumJS was used to create a pin located at the position of the Entity. The icons of GlobeLand30 were adopted for creating the pins. When a user clicks on a pin, the values of the significant properties of the JSON object that corresponds to the clicked pin are displayed in a window created using the InfoBox class of CesiumJS. The significant properties are the land cover class, the date the land cover data were submitted, the certainty of the user regarding the land cover class they stated and if available the comment of the user.
All the JSON objects that have the same land cover class are represented by an object created using the CustomDataSource class of CesiumJS so that they can be clustered. Clustering enables to get a less cluttered visualization, which enables to get insights regarding the distribution of land cover classes on the virtual globe for each zoom level. Each cluster is displayed with a pin created using the PinBuilder class of CesiumJS. The pin of each cluster of a land cover class has the colour of the pixels in the raster imagery that represents the clustered land cover class. Moreover, the number of Entities in each cluster is displayed on the pin. The clusters that are represented by pins on the virtual globe are displayed in Figure 13. Besides, this visualization enables to identify the differences between the VGI collected according to GlobeLand30 nomenclature and GlobeLand30 raster imagery visually. The differences suggest errors in GlobeLand30, which call upon further investigation. The data collected using the Land Cover Collector application can be used as reference data to perform the validation of GlobeLand30 computationally [79][80][81].

Mobility Visualization
Datasets of public transportation networks and traffic (i.e., number of vehicles) for each of 24 h of a day (10 October 2018) were generated by the team from Politecnico di Torino and stored on a GeoServer installation. They were retrieved using the WMS implementation of GeoServer and the ImageryLayer and WebMapServiceImageryProvider classes of CesiumJS were used to make the requests and place the retrieved images on the virtual globe. As traffic datasets are available for multiple time ranges, they can be visualized using the animation and timeline widgets of CesiumJS, using the method used in the applications related to ground deformation and LULC. However, it is worth to note that, in this application instead of ImageMosaic data store of GeoServer, 24 datasets for each city stored in vector format on GeoServer were used. In Figure 14, the visualization of traffic in Rome between 12:00 and 13:00 UTC is given.

Conclusions
Cartography has traditionally produced 2D maps. Globes, even though they provide the advantage of undistorted representation of the Earth, are uncomfortable to use. On the other hand, virtual globes are comfortable to use and allow to contextualize the represented phenomena more effectively compared to 2D maps. Furthermore, their potential is currently not fully explored. As a result, virtual globes were used for geovisualization in this work. For creating the virtual globes, open source software was used mainly so that they could be customized more freely compared to closed source alternatives.
Open standards were used to ensure the interoperability of the geoportal. Moreover, FOSS and open data were used for developing the geoportal primarily for diminishing the barriers to access it. The majority of the standards and FOSS used have been developed to deal with big geospatial data. As geospatial data have been steadily used more on the web and less with desktop applications, this research focused on the technological aspects for publication of geospatial data on the web using international and most commonly used open standards and the highest performing FOSS.
This work brings together the most recent web technologies and develops new software for multidimensional vector and raster geospatial data visualization and multidimensional raster geospatial data processing. As the challenges focused on in this research were examined not just theoretically, but also through developing applications, the applications can be considered emblematic for designing and implementing optimized solutions for tackling the aforementioned challenges. Furthermore, it is important to note that the application for processing through WCPS is only an example. The same software architecture can be applied to other operations that the WCPS standard allows executing.
This research can be continued by integrating geovisual analytics techniques, such as multiple-linked views into the geoportal. Moreover, the recently emerged 3D massive geospatial data streaming and rendering standard 3D Tiles can be used in the future if a FOSS alternative is available for creating datasets in this format. Furthermore, it would be beneficial to let the users write and execute their rasql queries through the geoportal. Additionally, the results of processing through WCPS that are in raster format can be visualized on the virtual globe. Furthermore, the users can be enabled to download all the processing results through the geoportal. Finally, the existing open data repositories can be explored, and the data that they host can be used in the geoportal without republishing them.