Commit 0d7acf50 authored by Luís de Sousa's avatar Luís de Sousa
Browse files

Completed tiling example with WCS

parent 2e126871
%% Cell type:markdown id: tags:
Tiling requests for large areas
=========================
The previous examples showed how to obtained an area of interest in a single request. In circumstances were the area of interest is large this method may not suceed. This can be the case with a large country like Brasil or a large geographic region like the Sahara. Requesting such a large map from the server can take too long and may even hit resource limits in the server.
The previous examples showed how to obtain an area of interest in a single request. In circumstances where the area of interest is large this method may not suceed. This can be the case with a large country like Brasil or a large geographic region like the Sahara. Requesting such a large map can take too long and may even hit resource limits in the server.
In these cases the best approach is to obtain the map with various requests. Each successive requests fetches a segment of the area of interest (a tile) than can later be composed together as a mosaic.
In these cases the best approach is to obtain the map with various requests. Each successive request fetches a segment of the area of interest (a tile) than can later be composed together as a mosaic.
Since the GetCoverage request will be issued several times, the first thing to do is to encapsulate it in a function:
Since the *GetCoverage* request will be issued several times, it is better to start by encapsulating it in a function:
%% Cell type:code id: tags:
``` python
from owslib.wcs import WebCoverageService
getTile(cov_id, subsets, resX, resY, format, tile):
def getTile(wcs, crs, cov_id, subsets, resX, resY, format, tile):
response = wcs.getCoverage(
identifier=[cov_id],
crs=crs,
subsets=subsets,
......@@ -37,10 +37,12 @@
wcs = WebCoverageService('https://maps.isric.org/mapserv?map=/map/phh2o.map',
version='2.0.1')
cov_id = 'phh2o_0-5cm_mean'
ph_0_5 = wcs.contents[cov_id]
crs = "http://www.opengis.net/def/crs/EPSG/0/152160"
xMin = -500000
xMax = 1500000
yMin = 2100000
yMax = 4200000
```
......@@ -52,11 +54,68 @@
%% Cell type:code id: tags:
``` python
resX = 250
resY = 250
tileSide = resX * 100
tileSide = resX * 2800
```
%% Cell type:markdown id: tags:
And the fun bit, a nested loop to obtain the tiles in succession. A counter, *tileN*, is used to name the tiles.
%% Cell type:code id: tags:
``` python
tileN = 0
subX = xMin
while subX < xMax:
subY = yMin
while subY < yMax:
subsets = [('X', subX, subX + tileSide),
('Y', subY, subY + tileSide)]
getTile(wcs, crs, cov_id, subsets, resX, resY,
ph_0_5.supportedFormats[0], str(tileN))
print("Saved tile " + str(tileN))
subY = subY + tileSide
tileN = tileN + 1
subX = subX + tileSide
```
%% Cell type:markdown id: tags:
[Index](index.ipynb) | [Previous](03-WCS-2.0.ipynb
It is more convinient to use the bundle of tiles as a single map. The VRT file format is the right tool for this job. You need to have [GDAL](https://gdal.org/) installed, and then use it directly from the command line:
%% Cell type:code id: tags:
``` python
!gdalbuildvrt ./data/Algeria_phh2o_0-5cm_mean.vrt ./data/phh2o_0-5cm_mean*.tif
```
%% Cell type:markdown id: tags:
The final map can now be plotted as in the previous example:
%% Cell type:code id: tags:
``` python
import rasterio
from rasterio import plot
ph = rasterio.open("./data/Algeria_phh2o_0-5cm_mean.vrt", driver="VRT")
%matplotlib inline
plot.show(ph, title='Mean pH between 0 and 5 cm deep in Algeria', cmap='gist_rainbow')
```
%% Cell type:markdown id: tags:
Experiment modifying the *tileSide*. Too small and it will result in a large number of tiles. Too large and it can overload the server.
%% Cell type:markdown id: tags:
[Index](index.ipynb) | [Previous](03-WCS-2.0.ipynb)
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment