5  Origin-Destination matrix

We created Origin and Destination (OD) points to create shortest path in each of the combinations using the costs from ORS instead of the distance.

For the centrality, we used the Global Human Settlement Built-up volume raster to obtain a weighted OD. Those regions with a higher built-up volume had a higher number of OD. The following map shows the weighted OD based on the building density that were used to calculate the core metropolitan centrality of the AoI (Core Metropolitan Area) and to calculate the intracity connectivity (Municipality).

library(sf)
library(DBI)

## Import Data from PostgreSQL database
municipalities_ghs <- st_read(connection,
        DBI::Id(schema="agile_gis_2025_rs",
                table="municipalities_ghs"))
build_masked <- terra::rast("~/agile-gscience-2024-rs-flood/data/source_data/build_masked.tif")
build_masked_na <- terra::subst(build_masked, from=0, to=NA)
core_metropolitan_area <- st_read(connection,
        DBI::Id(schema="agile_gis_2025_rs",
                table="core_metropolitan_area"))        
od_300_snapped_v2 <- st_read(
                              connection,
                              DBI::Id(schema="agile_gis_2025_rs",
                                      table="od_300_snapped_v2"))
## Create interactive map
### styles
library(raster)
library(leaflet)
pal <- colorNumeric(c("#faebdd", "#cd1c4e", "#110c24"), values(build_masked),
  na.color = "transparent")
pal_fct <- colorFactor(
  palette = c('#e41a1c', '#377eb8', '#4daf4a', '#984ea3', '#ff7f00','#ffff33','#a65628', '#f781bf','#999999' ),
  domain = od_200_each$shapename
)
### map
leaflet() |>
  addTiles()  |>
  addProviderTiles("CartoDB.Positron", group = "CartoDB.Positron")  |>
  addProviderTiles("Esri.WorldImagery", group = "Esri.WorldImagery")  |>                    addPolygons(data=municipalities_ghs,
                    fillColor = "yellow",
                    fillOpacity = 0.05,
                    weight = 1.5,
                    color = "red",
                    dashArray = "1",
                    group= "Municipalities") |> 
  addPolygons(data=core_metropolitan_area,
                    fillColor = "yellow",
                    fillOpacity = 0.05,
                    weight = 2,
                    color = "red",
                    dashArray = "1",
                    group= "Core Metropolitan Area") |> 
  addCircleMarkers(data=od_300_snapped_v2,
                    fillColor = "yellow",
                    fillOpacity = 0.8,
                    weight = 1.5,
                    radius= 3,
                    color = "black",
                    dashArray = "1",
                    group= "Weighted OD - C.Metropolitan Area") |> 
   addCircleMarkers(data=od_200_each,
                    color = ~pal_fct(shapename),
                    fillOpacity = 0.8,
                    weight = 1.5,
                    radius= 3,
                    dashArray = "1",
                    group= "Weighted OD - Municipality") |> 
  addRasterImage(build_masked_na,
                 opacity=0.7,
                 colors = pal,
                 group="Bulding density") |> 
     addLegend(pal = pal, values = values(build_masked),
    title = "Building Density (m³)") |> 
 addLayersControl(
    baseGroups = c("CartoDB.Positron", "Esri.WorldImagery"),
    overlayGroups = c("Municipalities","Core Metropolitan Area", "Weighted OD - C.Metropolitan Area","Weighted OD - Municipality","Bulding density"),
    options = layersControlOptions(collapsed = FALSE)
  ) |>
  hideGroup(c("Municipalities","Core Metropolitan Area","Weighted OD - Municipality")) |>
  setView( lng =-51.21566 , lat = -30.04176, zoom = 12 ) 

For the the redundancy, a component of the resiliene, we used a regular OD sampling in isochrones of 10 minutes for the hospitals equiped with ICU beds in the Area of Interest. The following map shows these Isochornes, hospitals as Point of Interest (PoI) and the regular OD samples.

## Import data:
regular_sampling_isochrones_snapped_ghs_aoi <- 
  sf::st_read(connection, Id(schema="agile_gis_2025_rs",table="regular_sampling_isochrones_snapped_ghs_aoi"))

## style:
factpal <- colorFactor(topo.colors(22), isochrones_poi_hospitals$ds_cnes)
## Map
leaflet() |>
  addTiles()  |>
  addProviderTiles("CartoDB.Positron", group = "CartoDB.Positron")  |>
  addProviderTiles("Esri.WorldImagery", group = "Esri.WorldImagery")  |>   
  addProviderTiles("CartoDB.DarkMatter", group = "CartoDB.DarkMatter") |> 
  addCircleMarkers(data=hospitals_poi,
                    fillColor = "yellow",
                    fillOpacity = 0.6,
                    weight = 3,
                    radius = 5,
                    color = "red",
                    dashArray = "1",
                    group= "Hospitals PoI") |> 
  addPolygons(data=isochrones_poi_hospitals,
                    fillOpacity = 0.2,
                    weight = 2,
                    color = ~factpal(ds_cnes),
                    dashArray = "1",
                    group= "Isochrones") |> 
  addCircleMarkers(data=regular_sampling_isochrones_snapped_ghs_aoi,
                    fillOpacity = 0.8,
                    weight = 1.5,
                    radius= 3,
                    color =~factpal(ds_cnes),
                    dashArray = "1",
                    group= "Regular OD") |> 
 addLayersControl(
    baseGroups = c("CartoDB.Positron", "Esri.WorldImagery","CartoDB.DarkMatter"),
    overlayGroups = c("Hospitals PoI","Isochrones", "Regular OD"),
    options = layersControlOptions(collapsed = FALSE)
  ) |> 
    setView( lng =-51.01134 , lat = -29.93259, zoom = 11.5 ) 

Troubleshootin 10 There were several tables when exporting od_300_v2 to GIS from PostgreSQL database. The reason is that unlike PostgreSQL, QGIS only admit one geometry type column. We solved this checking the attribute table in QGIS and choosing the table accordingly.

5.1 Weighted OD

5.1.1 Create weighted sampling points

5.1.1.1 For the Core Metropolitan Area

We masked the global dataset GHS-Built-up volume to the area of interest named Core Metropolitan Area and took 300 samples based the building density. A higher number of points caused the ERROR: invalid memory alloc request size 2450771848, however, a possible solution is to iterate using batch queries. Although the limitation of the OD to 300 points could be an underrepresented of the phenomena, we did observed similar patterns in the centrality with visual inspections. The function spatSample from the library terra (Hijmans 2024) created the table od_300_v2 that is used to calculate the core metropolitan centrality.

# Import data
library(sf) # handle vector files
library(terra) # handle raster files
buildings <- terra::rast("/home/ricardors/agile-gscience-2024-rs-flood/data/source_data/GHS_BUILT_V_E2020_GLOBE_R2023A_4326_100_V1_0_RioGrandeDoSul.tif")
pop_ghs <- 300 # set up number of sampling points
## Select area of interest from the database
municipalities_ghs <- sf::read_sf(connection,
                                  DBI::Id(schema="agile_gis_2025_rs",
                                          table="municipalities_ghs"))
municipalities_ghs_reproj <- municipalities_ghs |>
                                            st_transform(crs(buildings))
build_cropped <- crop(buildings,
                      municipalities_ghs_reproj)
## Reproject for sampling
build_4326 <- terra::project(build_cropped,
                             crs(municipalities_ghs))
build_masked <- terra::mask(build_4326,
                            municipalities_ghs)
## Weighted sampling
set.seed(123)
od <- spatSample(build_masked,
                     pop_ghs,
                    "weights",
                     as.points = TRUE) |>
                st_as_sf() |>
                st_transform(4326)
## Rename column
names(od)[1] <- "building"
## Upload the sampling points to the database
sf::dbWriteTable(connection, DBI::Id(schema = "agile_gis_2025_rs",
table = "od_300_v2"), od)

# time: 721 ms

building	id	geometry
47979	1	POINT (-50.99947 -29.92425)
26943	2	POINT (-51.14501 -29.86022)
41710	3	POINT (-51.11673 -30.0149)
18807	4	POINT (-51.13919 -30.0174)
9664	5	POINT (-51.06766 -29.94671)

5.1.1.2 For the Intracity connectivity

We iterated through each municipality masking the GHS-Built-up volume and taking 200 weighted OD samples. The list that stored all the iterations is converted into a spatial dataframe (sf).

# Load the municipalities 
municipalities_ghs <- sf::st_read(connection,
                                  Id(schema="agile_gis_2025_rs",
                                     table="municipalities_ghs"))
# Select size of sample population
n_samples <- 200

# Create a list before looping for each municipality to generate sampling points
sampled_points_list_200 <- list()
library(tidyverse)
for (idx in 1:nrow(municipalities_ghs)){
  shapename <- municipalities_ghs$shapename[idx]
  municipality_geom <- municipalities_ghs$wkb_geometry[idx]
  n_od_samples <- n_samples
  masked_ghs_aoi <- mask(build_masked, vect(municipality_geom))
  
od <- spatSample(masked_ghs_aoi, size = n_od_samples, method = "weights", as.points = TRUE) 
 od_sf <- st_as_sf(od) |> st_transform(4326)
  
  # df with the point id and muncipality (shapename)
  od_sf <- od_sf |> 
    mutate(
      id_sample = row_number(),
      shapename = shapename
    ) |> 
    select(id_sample, shapename, geometry = geometry)
  
  # list that stores all the points
  sampled_points_list_200[[idx]] <- od_sf
}
# Create a dataframe from the previous list created during the loop
sampled_points_df_200 <- bind_rows(sampled_points_list_200)
# Export from R to PostgreSQL
DBI::dbWriteTable(connection, Id(schema="agile_gis_2025_rs","od_200_each"),sampled_points_df_200)
# time: 1468 ms

The following table shows the first 5 observations from the table od_200_each where the municipality is indicated in the field shapename.

id_sample	shapename	geometry
1	Alvorada	POINT (-51.08014 -29.98829)
2	Alvorada	POINT (-51.0685 -30.02239)
3	Alvorada	POINT (-51.08097 -30.03902)
4	Alvorada	POINT (-51.05436 -30.0149)
5	Alvorada	POINT (-51.08097 -29.99577)

5.1.2 Snap weighted sampling points

5.1.2.1 For the Core Metropolitan Area

The sampling points of the OD matrix are not necessary located on the road network, therefore it is necessary to snap these points. Previous to snap the points, we made the largest component of the road network routable. (i.e. ghs_ao_net_largest). Subsequently, we indexed the geometries and ID variables to improve computational efficiency. We calculated the Euclidean distance between each sampling point and the nodes from the road network nodes. By limiting the results by one, we selected the closest node for each sampling point.

--- (optional) Create "id" column as a pre-requisite for snapping if it is not created before.
ALTER TABLE od_300_v2
ADD COLUMN id INTEGER PRIMARY KEY GENERATED ALWAYS AS IDENTITY;

--- Similar to [ORN-3], we make the largest component routable
---   First we cast the variables according to pgrouting requirements
ALTER TABLE ghs_aoi_net_largest 
    ALTER COLUMN "toId" type bigint,
    ALTER COLUMN "fromId" type bigint,
    ALTER COLUMN "id" type bigint;
--- Second we create the auxiliary table ghs_aoi_net_largest_vertices_pgr
SELECT pgr_createverticesTable('ghs_aoi_net_largest',
the_geom:='geom',
source:='fromId',
target:='toId');

--- Create  indexes to improve queries performance.
CREATE INDEX idx_od_300_v2_id ON od_300_v2 USING btree(id);
CREATE INDEX idx_od_300_v2_geom ON od_300_v2 USING gist(geometry);
CREATE INDEX idx_ghs_aoi_net_largest_vertices_pgr ON ghs_aoi_net_largest_vertices_pgr USING btree(fid);
--- Snap the OD to the nodes of ghs_aoi_net_largest_vertices_pgr
EXPLAIN ANALYZE
CREATE TABLE od_300_snapped_v2 AS
SELECT DISTINCT ON (net.id)
  pt.id AS pt_id, --- original position of the Origin-Destination
  net.id AS snapped_id, --- snapped position of the Origin-Destination
  net.the_geom
FROM
(SELECT * FROM od_300_v2 AS pt) as pt
CROSS JOIN
LATERAL (SELECT
            *
        FROM
          ghs_aoi_net_largest_vertices_pgr AS net
        ORDER BY net.the_geom <->pt.geometry
        LIMIT 1) AS net;
---time: 35.474 ms

5.1.2.2 For the Intracity connectivity

The OD matrix is limited to each municipality instead of the core metropolitan area. We created a numeric column named shapename_id to iterate each municipalities’ networks.

--- Create a ID for each municipality
CREATE TABLE od_200_each_group_id AS
SELECT 
    *,
dense_rank() OVER (ORDER BY shapename) 
FROM od_200_each;
---time: 13.858 ms

First, we create an empty table that contained the final columns.

library(glue)
library(DBI)
od_200_each_group_id <- st_read(connection,"od_200_each_group_id")

table_snapped <- glue_sql("CREATE TABLE table_snapped 
                (pt_id numeric,
                snapped_id numeric,
                shapename_id varchar(200),
                shapename varchar(200),
                geom geometry)", .con =connection)
table_snapped_db <- dbGetQuery(conn=connection, table_snapped )

Then, we snapped the sampling points to the closest vertice of the network. The table net_largest_1_vertices_pgr was created in the loop shown in code chunk [ORN-6]. The final table table_snapped contained the sampling point weighted on the building density for each of the municipalities.

tic()  
for(idx in 1:9){
  tic(glue("run {idx}/9")) 
  table_name_1_centrality <- glue("net_largest_{as.character(idx)}")
  temp_table_name <- glue("temp_centrality_{idx}")
  table_bidirect_id <- glue("t_bidirect_id_{idx}")
  dbExecute(glue_sql("DROP TABLE IF EXISTS {`temp_table_name`}", .con=connection),
          conn=connection)  
query_1_central <- glue_sql("CREATE TABLE {`temp_table_name`} AS
 SELECT   
 b.id,
 b.geom,
 count(b.geom) AS centrality 
 FROM  pgr_dijkstra('SELECT id,
                            fromid AS source,
                            toid AS target,
                            weight AS cost
                      FROM {`table_name_1_centrality`}',
                      ARRAY(SELECT 
                        snapped_id AS start_id 
                        FROM table_snapped 
                        WHERE shapename_id = {idx}::text),
                      ARRAY(SELECT
                        snapped_id AS end_id 
                        FROM table_snapped
                        WHERE shapename_id = {idx}::text),
                      directed := TRUE) j
                      LEFT JOIN {`table_name_1_centrality`} AS b
                      ON j.edge = b.id
                      GROUP BY  b.id, geom;", .con=connection)

dbExecute(con=connection, query_1_central)
toc()
}
toc()

# Time: 53805 ms

Troubleshooting 9 The query did not find the required tables. Since this iteration required the tables net_largest_x, where x is one of the municipalities id, we observed that these tables were in the public schema. We solved this by defining the search_path or the schema.

5.2 Regular OD

The subject of the study that uses a regular OD is the hospital accessibility and their lack of redundancy. First we create the field location using isochrones generated with ORS. Then we create regular samples and finally we snap these OD to the road networks’ nodes.

5.2.1 Select Point of Interest (PoI)

We used ORS to create 10 minutes ischrones for each of our point of interest, the hospitals. Before selecting running the PostgreSQL query, we created indexes in the geometry column and id column. The first common table expressin (CTE) all_hospitals included the non-spatial information from the hospitals with their location joining both tables by the hospital code cne. Subsequently, the CTE poi_hospitals filtered all the hospitals selecting only those that intersected with the core metropolitan area. Lastly, the table ghs_aoi_net_largest_vertices_pgr created by pgrouting provided the nodes or vertices to which the hospitals are snapped.

CREATE INDEX idx_etlcnes_hospital_selected_cnes ON etlcnes_hospital_selected USING btree("CNES");
CREATE INDEX idx_hospitals_cd_cnes ON hospitals USING btree(cd_cnes);
EXPLAIN ANALYZE
CREATE TABLE hospitals_poi AS
--- Add ICU bed capacity from ETLCNES to the hospitals
WITH all_hospitals AS (SELECT DISTINCT ON (cd_cnes)
    cd_cnes,
    ds_cnes,
    wkb_geometry AS geom,
        "QTLEITP1",
    "QTLEITP2",
    "QTLEITP3",
    "NAT_JUR",
    tp_risco,
    nm_razao_s,
    nm_municip
FROM
    hospitals AS geom_hospital
LEFT JOIN
    etlcnes_hospital_selected AS df_hospital
ON
    df_hospital."CNES"=geom_hospital.cd_cnes),
--- Define Point Of Interest (POI) using Area of Interest (AoI)
     poi_hospitals AS (
SELECT 
    h.*
FROM 
    all_hospitals AS h,
    core_metropolitan_area AS poi
WHERE st_intersects(h.geom, poi.geom))
--- Add the original geometry and the snapped position
SELECT DISTINCT ON (h.cd_cnes)
    h.*,
    f.the_geom <-> h.geom AS distance,
    h.geom AS geom_hospital,
    f.the_geom AS geom_node
FROM poi_hospitals h
LEFT JOIN LATERAL
(SELECT 
    id, 
    the_geom
FROM ghs_aoi_net_largest_vertices_pgr AS net
ORDER BY
    net.the_geom <-> h.geom
LIMIT 1) AS f ON true

--- Time: 814.451 ms

5.2.2 Create Isochrones from PoI

The ORS API¹ created the area that is reached in 10 minutes from each of the hospitals, where we took the regular OD.

Error: Openrouteservice API request failed [3004] Parameter ‘locations=22’ is out of range. Maximum possible value is 5. The function ors_api_key registered a personal token required to run the ORS services by the API. Then, the arguments of the function ors_isochrones took the first 5 hospitals location, the range and interval of 600 seconds (10 minutes) and the output format, a sf object. The reason to repeat this process five times was that the maximum possible values for location was 5. This approach is quite rudimentary, yet it worked. The last step was to remove the “center” column that resulted on an non-valid format to copy the R object into PostgreSQL using the function st_write.

library(sf)
library(openrouteservice)
# Use the free API from ORS replacing "<>" for the key value.
openrouteservice::ors_api_key("<https://openrouteservice.org/dev/#/login?redirect=https://openrouteservice.org/dev/%23/api-docs>")

hospitals_poi <- sf::st_read(connection, DBI::Id(schema="agile_gis_2025_rs" ,"hospitals_poi"))
library(tidyverse)
## Part I: PoI (1-5) 
tic()
isochrone_part_1 <- 
  ors_isochrones(st_coordinates(hospitals_poi$geom_node[1:5]),
                                range = 600,
                                interval = 600,
                                output = "sf") |>
                            mutate(ds_cnes = hospitals_poi$ds_cnes[1:5],
                                   cd_cnes = hospitals_poi$cd_cnes[1:5],
                                   id = hospitals_poi$id[1:5])
## Part II: PoI (6-10)
isochrone_part_2 <- 
  ors_isochrones(st_coordinates(hospitals_poi$geom_node[6:10]),
                                range = 600,
                                interval = 600,
                                output = "sf") |>
                            mutate(ds_cnes = hospitals_poi$ds_cnes[6:10],
                                   cd_cnes = hospitals_poi$cd_cnes[6:10],
                                   id = hospitals_poi$id[6:10])
## Part III: PoI (11-15)
isochrone_part_3 <- 
  ors_isochrones(st_coordinates(hospitals_poi$geom_node[11:15]),
                                range = 600,
                                interval = 600,
                                output = "sf") |>
                            mutate(ds_cnes = hospitals_poi$ds_cnes[11:15],
                                   cd_cnes = hospitals_poi$cd_cnes[11:15],
                                   id = hospitals_poi$id[11:15])
## Part IV: PoI (16-20)
isochrone_part_4 <- 
  ors_isochrones(st_coordinates(hospitals_poi$geom_node[16:20]),
                                range = 600,
                                interval = 600,
                                output = "sf") |>
                            mutate(ds_cnes = hospitals_poi$ds_cnes[16:20],
                                   cd_cnes = hospitals_poi$cd_cnes[16:20],
                                   id = hospitals_poi$id[16:20])
## Part V: PoI (20-22)
isochrone_part_5 <- 
  ors_isochrones(st_coordinates(hospitals_poi$geom_node[21:22]),
                                range = 600,
                                interval = 600,
                                output = "sf") |>
                            mutate(ds_cnes = hospitals_poi$ds_cnes[21:22],
                                   cd_cnes = hospitals_poi$cd_cnes[21:22],
                                   id = hospitals_poi$id[21:22])
## Merge 5 parts

isochrones_poi_hospitals <- rbind(isochrone_part_1,
                                  isochrone_part_2,
                                  isochrone_part_3,
                                  isochrone_part_4,
                                  isochrone_part_5 )  |>
                                select(-center) 
# the variable ""center"" is a list. Writting will return  no applicable method for st_write.
toc()
# time: 1925 ms

5.2.3 Create regular sampling points

We used a function in PostGIS to sample regularly² the OD to assesss the hospitals accessibility and lack of redundancy. In this case, the function is limited to the Spatial Reference IDentifier (SRID) 4326. Since this SRID use degrees and their distance depends on the distance from the equator, the accuracy could be improved by using another SRID adapted to the project. We used this SRID because the accuracy between the OD were not as important as using a SRID that could be used on a global scale.

CREATE OR REPLACE FUNCTION I_Grid_Point_Series(geom geometry, x_side decimal, y_side decimal, spheroid boolean default false)
RETURNS SETOF geometry AS $BODY$
DECLARE
x_max decimal;
y_max decimal;
x_min decimal;
y_min decimal;
srid integer := 4326;
input_srid integer;
x_series DECIMAL;
y_series DECIMAL;
BEGIN
CASE st_srid(geom) WHEN 0 THEN
  geom := ST_SetSRID(geom, srid);
  RAISE NOTICE 'SRID Not Found.';
    ELSE
        RAISE NOTICE 'SRID Found.';
    END CASE;

    CASE spheroid WHEN false THEN
        RAISE NOTICE 'Spheroid False';
    else
        srid := 900913;
        RAISE NOTICE 'Spheroid True';
    END CASE;
    input_srid:=st_srid(geom);
    geom := st_transform(geom, srid);
    x_max := ST_XMax(geom);
    y_max := ST_YMax(geom);
    x_min := ST_XMin(geom);
    y_min := ST_YMin(geom);
    x_series := CEIL ( @( x_max - x_min ) / x_side);
    y_series := CEIL ( @( y_max - y_min ) / y_side );
RETURN QUERY
SELECT st_collect(st_setsrid(ST_MakePoint(x * x_side + x_min, y*y_side + y_min), srid)) FROM
generate_series(0, x_series) as x,
generate_series(0, y_series) as y
WHERE st_intersects(st_setsrid(ST_MakePoint(x*x_side + x_min, y*y_side + y_min), srid), geom);
END;
$BODY$ LANGUAGE plpgsql IMMUTABLE STRICT;

The following code created a bounding box required to run the created function in the previous step.

CREATE TABLE core_metropolitan_bbox AS
SELECT 
    ST_Extent(geom)::geometry AS geom
FROM 
    core_metropolitan_area;

Based on the isochrones of the healthcare facilities created from ORS, we regularly sampled the OD separated by 0.00912º.

EXPLAIN ANALYZE
CREATE TABLE  regular_sampling_isochrones_ghs_aoi AS
WITH multipoint_regular AS(
SELECT 
    I_Grid_Point_Series(geometry, 0.00912,0.00912, false) AS geom,
    ds_cnes,
    cd_cnes
FROM 
    isochrones_poi_hospitals),
point_regular AS(
SELECT 
    st_setsrid((st_dump(geom)).geom, 4326)::geometry(Point, 4326) AS geom,
    ds_cnes,
    cd_cnes
FROM  multipoint_regular)
SELECT 
    cd_cnes,
    ds_cnes,
    geom,
    row_number() over() AS id
FROM 
    point_regular; 
--- Time: 6083.595 ms

5.2.4 Snap regular sampling points

The code chunk ACR-9 snapped the OD to the closest node of the road network of each hospital isochrone. We noticed that the final number of OD may be slightly different from the rmarkdown_data of the github repository. This difference was not determinant to draw different conclusions from ours, therefore we acknowledge them and do not consider them.

EXPLAIN ANALYZE
CREATE TABLE regular_sampling_isochrones_snapped_ghs_aoi AS
SELECT DISTINCT ON (net.id)
    net.id,
    ds_cnes,
    cd_cnes,
    net.the_geom  AS geom
FROM regular_sampling_isochrones_ghs_aoi AS pt
CROSS JOIN LATERAL 
    (SELECT 
        *
    FROM  ghs_aoi_net_largest_vertices_pgr AS net
    ORDER BY
    net.the_geom <-> pt.geom
    LIMIT 1) AS net;        
--- Time: 166.008 ms

Hijmans, Robert J. 2024. Terra: Spatial Data Analysis. https://rspatial.org/.

---

1. https://openrouteservice.org/dev/#/login?redirect=https://openrouteservice.org/dev/%23/api-docs

2. https://gis.stackexchange.com/questions/4663/creating-regular-point-grid-inside-polygon-in-postgis