using GMT, RemoteS

# Tell GMT to display the images in VSCode
ENV["DISPLAY_IN_VSC"] = "a";

Definições:

• Um tile é sempre uma imagem com 256 x 256 pixeis
• Cada tile tem dois nomes únicos. Um baseado num esquema que se chama uma àrvore de quadrados e constituido apenas pelos algarismos 0,1,2 e 3. O número total de algarismos do nome indica o seu nível de zoom. No outro esquema, chamado XYZ, o nome é construido a partir do número da linha e da coluna de uma certa matriz (mais abixo veremos exemplos).
• A área coberta por um tile depende do seu nível de zoom. Por exemplo, ao nível 21, um tile cobre um rectângulo com o lado vertical com ~15 m.
• gratícula, ou se quisermos, malha, são as linhas horizontais e verticais (a vermelho nas imagens seguintes) que delimitam cada tile.
• Uma imagem devolvida pela função mosaic é constituida por um ou uma junção de vários tiles colados lado-a-lado.

Uma imagem com um zoom de nível 14 cobre uma área rectangular com perto de 2 km de lado N-S. Vejamos um exemplo da zona das Gambelas.

I14 = mosaic(-7.966902, 37.042682, provider=:goog, zoom=14);

viz(I14)

Esta imagem tem uma resolução muito baixa e apenas servirá para nos guiar. Um aspecto muito interessante destas imagens (imagens dos TileMapServers) é que cada tile tem dois nomes únicos que podemos nos atributos ver se usarmos a opção mesh=true

D14 = mosaic(-7.966902, 37.042682, zoom=14, mesh=true)

Comment:	["Zoom level = 15"]
Attribute table
┌────────────────┬──────────────┐
│       quadtree │          XYZ │
├────────────────┼──────────────┤
│ 03311032210321 │ 7829,6374,14 │
└────────────────┴──────────────┘
BoundingBox: [-7.97607421875, -7.9541015625, 37.03763967977139, 37.05517710666082]
PROJ: +proj=longlat +datum=WGS84 +units=m +no_defs
Header:	 -W

5×2 GMTdataset{Float64, 2}
 Row │      Lon      Lat 
─────┼───────────────────
   1 │ -7.97607  37.0376
   2 │ -7.97607  37.0552
   3 │ -7.9541   37.0552
   4 │ -7.9541   37.0376
   5 │ -7.97607  37.0376

Os nomes deste tile são "03311032210321" e "7829,6374,14" que em si têm codificado as coordenadas da imagen. No que se segue, vamos usar a codificação XYZ. O 14 no nome representa o nível de zoom (o X e o Y são as linhas e colunas de uma certa matriz).

Podemos então usar estes nomes para pedir mais imagens (ou mais gratículas) à função mosaic. Vamos começar por pedir as gratículas do nível 15 e sobrepô-las à imagem I14.

D15 = mosaic("7829,6374,14", zoom=1, mesh=true)

Vector{GMTdataset} with 4 segments
Showing first segment. To see other segments just type its element number. e.g. D[2]

Comment:	["Zoom level = 16"]
Attribute table
┌─────┬─────────────────┬────────────────┐
│ Row │ quadtree        │ XYZ            │
├─────┼─────────────────┼────────────────┤
│   1 │ 033110322103210 │ 15658,12748,15 │
│   2 │ 033110322103212 │ 15658,12749,15 │
│   3 │ 033110322103211 │ 15659,12748,15 │
│   4 │ 033110322103213 │ 15659,12749,15 │
└─────┴─────────────────┴────────────────┘
BoundingBox: [-7.97607421875, -7.965087890625, 37.046408899699564, 37.05517710666082]
Global BoundingBox: [-7.97607421875, -7.9541015625, 37.03763967977139, 37.05517710666082]
PROJ: +proj=longlat +datum=WGS84 +units=m +no_defs
Header:	 -W

5×2 GMTdataset{Float64, 2}
 Row │      Lon      Lat 
─────┼───────────────────
   1 │ -7.97607  37.0464
   2 │ -7.97607  37.0552
   3 │ -7.96509  37.0552
   4 │ -7.96509  37.0464
   5 │ -7.97607  37.0464[-7.97607421875, -7.965087890625, 37.046408899699564, 37.05517710666082]
Global BoundingBox: [-7.97607421875, -7.9541015625, 37.03763967977139, 37.05517710666082]
PROJ: +proj=longlat +datum=WGS84 +units=m +no_defs
Header:	 -W

5×2 GMTdataset{Float64, 2}
 Row │      Lon      Lat 
─────┼───────────────────
   1 │ -7.97607  37.0464
   2 │ -7.97607  37.0552
   3 │ -7.96509  37.0552
   4 │ -7.96509  37.0464
   5 │ -7.97607  37.0464

Podemos ver que agora temos 4 quadrados pois aumentar um nível de zoom divide a área por dois e como pedimos os mesmos limites temos 4 vezes mais quadrados. Notem bem a novidade the agora termos usado "7829,6374,14", zoom=1. Ou seja, como delimitador da região usámos o mome to tile I14, e como zoom usámos 1, que significa adicionar 1 ao nível determinado pelo uso do nome 7829,6374,14.

Vejamos o que se obtem sobrepondo a gratícula e a image I14. E vamos usar mais uma novidade que é a de dizer que queremos os nomes dos tiles impressos na figura.

viz(I14, plot=(data=D15, lc=:red, lw=0.5, labels="attrib=XYZ"))

E agora podemos aplicar este princípio para obter uma imagem exactamente da mesma zona mas com un nível de zoom superior. Por exemplo para obter uma imagem de nível 15, fazemos:

I15 = mosaic("7829,6374,14", provider=:google, zoom=1);
viz(I15, plot=(data=D15, lc=:red, lw=0.5, labels="attrib=XYZ"))

Estão a ver que para chegarmos ao nível 19 é só continuar o processo, mas com cuidado porque a cada nível a mais quadruplicamos o número de tiles. Como exemplo vamos ver o que dá a gratícula do nível 17 sobreposta na imagem I15

D17 = mosaic("7829,6374,14", zoom=3, mesh=true);
viz(I15, plot=(data=D17, lc=:red, lw=0.5, labels="attrib=XYZ"))

Ou seja, para se chegar ao nível 19 basta-nos dar o nome de um tile 17 e somar 2 ao zoom. Por exemplo, para obtermos uma imagem semelhante à que usámos nas aulas vamos usar o tile "62634,50997,17" como base.

I17_19 = mosaic("62634,50997,17", zoom=2, provider=:goog);
viz(I17_19)

Notem que temos de ler o nome to tile na figura e depois escrevê-lo à mão, o que é sempre uma possível fonte de erro. A alternativa é usar a tabela das gratículas sabendo que elas estão numeradas por ao longo das colunas. Quer isto dizer que tendo o D17 64 quadrados e o que escolhemos foi o quadrado 22 (começa a contar no quadrado do canto superior esquerdo e vai por aí abaixo).

D17[22].attrib

Dict{String, Union{String, Vector{String}}} with 2 entries:
  "quadtree" => "03311032210321212"
  "XYZ"      => "62634,50997,17"

onde se confirma que o tile usado foi o "62634,50997,17"

Claro que já estou a ouvir os lamentos: ena tantos tiles para fazer. Sim, talvez mas muitos vêm-se logo que não têm piscinas e depois este foi um exemplo para fazer imagens de nível 19 que sejam pequenas e possíveis de manipular em computadores menos artilhados. Um númeo de imagens de nível 19 mais reduzido poderá ser obtido se se usar como tile de base os de nível 16, com os quais se pedem mosaicos 3 níveis mais acima. E quem se quiser atrever até pode experimentar a fazer a identificação das piscinas ao nível 20 em vez do 19.

Um problema, contudo, com exemplo apresentado na imagem acima é que ela possivelmente tem demasiado poucas piscinas para nela podermos desenhar (no Mirone, ou QGis, ou outro) os poligonos das classes que vamos usar como treino. Por isso vamos ver como obter uma área maior mas com o zoom de 19.

D16 = mosaic("7829,6374,14", zoom=2, mesh=true);
viz(I15, plot=(data=D16, lc=:red, lw=0.5, labels="attrib=XYZ"))

E agora a image maior da área das Gambelas com um zoom de 19 obtem-se somando 3 ao tile "31317,25498,16"

I16_19 = mosaic("31317,25498,16", zoom=3, provider=:goog);
viz(I)

O que eu me referia na questão das dimensões das imagens e como isso pode ser pesado para alguns computadores é que as imagens _I17_19_ e _I16_19_ têm dimensões:

size(I17_19)

(1139, 909, 3)

size(I16_19)

(2278, 1818, 3)

Não é que a última seja particularmente grande mas internamente no processo de treino/classificação há cópias que precisam de ser feitas. Se poderem usar o esquema das imagens maiores terão menos trabalho.