Rasterio: Numpy for Geospatial Data

Olá, daqui é o Alex da DEV STORIES DOT EU. Rasterio: Numpy para Dados Geoespaciais, episódio 1 de 5. Queres extrair valores de pixels de uma imagem de satélite, mas as libraries de imagem standard engasgam-se com formatos geoespaciais, e os bindings oficiais dão a sensação de estar a escrever código C em Python. Essa tensão é exatamente o motivo pelo qual o Rasterio existe. O Rasterio é uma library que lê e escreve dados geoespaciais baseados em grelha, como ficheiros GeoTIFF. É muito provável que já tenhas ouvido falar do GDAL, a library geoespacial standard da indústria. Um equívoco comum é achar que o Rasterio substitui o GDAL. Não substitui. Na verdade, o Rasterio usa o GDAL under the hood. A diferença está na interface. Se alguma vez usaste os bindings nativos de Python do GDAL, sabes que eles expõem uma API C desajeitada, cheia de gestão manual de recursos e uma sintaxe que parece totalmente estranha ao Python. O Rasterio deixa isso para trás. Ele fornece uma interface moderna e Python-first. Ele trata os datasets raster como ficheiros normais e lida com os seus pixels como arrays Numpy standard. Supõe que tens um GeoTIFF de imagens Landsat no teu disco. Abres este ficheiro usando a função rasterio ponto open. Como o Rasterio adota os patterns standard de Python, fazes isso usando um with statement. Isto cria um context manager. Tal como ao abrir um ficheiro de texto básico, o bloco with garante que o dataset é fechado de forma limpa e a memória é libertada no momento em que o teu código termina de correr. Nunca tens de destruir objetos manualmente, nem de te preocupar com memory leaks de ficheiros mal fechados. Dentro desse bloco with, a função open dá-te um objeto dataset reader. Antes de leres quaisquer pixels reais, podes inspecionar o ficheiro. Podes pedir ao dataset o seu count, que te diz o número total de bandas, ou layers, na imagem. Para uma imagem Landsat, este count pode ser algo como sete ou onze, representando diferentes comprimentos de onda de luz. Podes ler as propriedades width e height para obteres as dimensões espaciais em pixels. Também podes verificar a propriedade dtypes. Isto dá-te o data type dos pixels para cada banda, como unsigned 16-bit integers ou 32-bit floats. Toda esta metadata está disponível instantaneamente, sem carregar os dados pesados da imagem para a memória do teu sistema. Aqui está o ponto chave. Quando estiveres pronto para olhar para os pixels reais, chamas o método read no objeto do dataset. Podes passar um index específico para carregar apenas uma layer. Presta atenção a este detalhe: as bandas do Rasterio são one-indexed, o que significa que pedes a banda um, e não a banda zero. Se chamares o read e passares o número um, o método devolve esses pixels como um array Numpy bidimensional standard. Se chamares o read sem quaisquer argumentos, ele lê tudo, devolvendo um array Numpy tridimensional que contém todas as bandas. Não há nenhum objeto pixel especial ou proprietário. Obténs simplesmente um array de números. Assim que tiveres esse array Numpy, estás de volta a território familiar. Podes fazer slice, correr operações estatísticas, ou alimentá-lo diretamente em modelos de machine learning usando exatamente o mesmo código que usas para qualquer outra matriz. O Rasterio atua como uma ponte. Ele lida com os formatos de ficheiro e metadados complicados no disco, e entrega-te um objeto matemático limpo em memória. Se gostas do programa e queres apoiar o que fazemos, podes procurar por DevStoriesEU no Patreon — é uma ajuda enorme. O verdadeiro poder do Rasterio não é reinventar o processamento geoespacial, mas sim torná-lo aborrecido; ao transformar imagens de satélite complexas em arrays Numpy standard, permite-te parar de lutar contra formatos GIS e começar a fazer data science a sério. Obrigado por ouvirem, e happy coding para todos!

Olá, daqui fala o Alex da DEV STORIES DOT EU. Rasterio: Numpy para Dados Geoespaciais, episódio dois de cinco. Um array Numpy básico é apenas uma grelha plana de números. Se o visualizares, obténs simplesmente um retângulo colorido no teu monitor. O que transforma realmente esses números abstratos num mapa preciso do Planalto do Colorado? A resposta é a Transformação Afim. As pessoas costumam achar que os dados geoespaciais são apenas uma imagem onde os cantos estão vagamente fixados às coordenadas do mapa. Não é bem assim. A transformação é uma ligação matemática rigorosa que calcula a posição exata no mundo real de cada pixel do teu array. Antes de poderes mapear pixels na Terra, tens de saber como estás a medir a Terra. Este é o Coordinate Reference System. Lês isto do dataset acedendo a dataset ponto crs. Esta propriedade diz ao Rasterio se as tuas coordenadas espaciais são medidas em graus de longitude e latitude, ou em metros em relação a um equador e a um meridiano principal específicos. Assim que o sistema de referência estiver definido, precisas da lógica de mapeamento. Encontras isto verificando dataset ponto transform. Isto devolve uma matriz de transformação afim. Isto soa a álgebra linear pesada, mas na prática é um conjunto de seis números que descrevem a área física da tua grelha. A matriz define as coordenadas espaciais exatas x e y do canto superior esquerdo do pixel do topo esquerdo. Define a largura física de um pixel, por exemplo, exatamente trinta metros. Define a altura física de um pixel, que geralmente é negativa porque as linhas da imagem contam para baixo enquanto as coordenadas do mapa contam para cima. A matriz também guarda valores de rotação para o caso de o satélite estar inclinado quando a imagem foi capturada. Aqui está o ponto chave. Nunca tens de multiplicar estes valores da matriz à mão. O Rasterio dá-te métodos diretos para saltar para trás e para a frente entre o espaço de pixels e o espaço de coordenadas. Supõe que estás a analisar uma imagem e encontras uma anomalia num local específico, digamos, na linha quinhentos e coluna mil. Para descobrires onde isso fica no planeta real, passas esses dois números inteiros para o método dataset ponto xy. O Rasterio passa a linha e a coluna pela matriz afim e devolve as coordenadas espaciais exatas x e y do ponto central desse pixel específico. Também podes correr esta máquina ao contrário. Digamos que tens as coordenadas espaciais exatas de uma localização física. Talvez seja um desfiladeiro localizado exatamente a cem quilómetros a leste e cinquenta quilómetros a sul da origem da tua imagem. Precisas de extrair os dados de pixel subjacentes para esse desfiladeiro. Pegas nos teus valores espaciais x e y e passas-os para o método dataset ponto index. O Rasterio inverte a matriz afim, processa as tuas coordenadas e devolve os números inteiros exatos da linha e da coluna. Depois usas esses números inteiros para fazer um slice no teu array Numpy e ler a medição exata registada pelo satélite. A transformação afim isola o teu array de dados brutos da tua geometria espacial, garantindo que, não importa como cortes ou leias os teus dados, nunca perdes a tua posição exata no globo. É tudo por hoje. Obrigado por ouvires — vai construir algo fixe.

Olá, daqui é o Alex da DEV STORIES DOT EU. Rasterio: Numpy para Dados Geoespaciais, episódio 3 de 5. Se tentares ler um GeoTIFF de dez gigabytes diretamente para a memória do teu portátil, o teu script vai quase de certeza crashar. Podes precisar apenas de um pequeno pedaço de pixels para testar um algoritmo, mas o sistema tenta engolir o continente inteiro de uma só vez. A solução para este limite de memória é o Windowed Processing. O Windowed Processing permite-te ler ou escrever subsets retangulares específicos de um ficheiro raster, enquanto o resto do ficheiro fica guardado em segurança no teu disco rígido. Isto é totalmente baseado em índices. Hoje não estamos a lidar com coordenadas geográficas. Estamos a lidar estritamente com linhas e colunas de pixels. Para leres um subset, precisas de uma forma de dizer ao Rasterio que pixels capturar. Fazes isto ao importar o objeto Window do módulo windows do rasterio. A forma standard de definir uma Window é fornecendo quatro números: column offset, row offset, width e height. Presta atenção a esta parte. A ordem dos argumentos importa. Exige primeiro o column offset e, a seguir, o row offset. O column offset é a tua posição horizontal inicial, medida a partir da margem esquerda. O row offset é a tua posição vertical inicial, medida para baixo a partir da margem superior. Depois dos offsets, especificas quantos pixels de largura e quantos pixels de altura a tua seleção deve ter. Considera um ficheiro de teste enorme, com o tamanho de um gigabyte. Queres extrair um pequeno tile com 256 pixels de largura e 512 pixels de altura. Queres que este tile comece a mil pixels da esquerda e a dois mil pixels do topo. Instancias a tua Window passando mil para o column offset, dois mil para o row offset, 256 para a width e 512 para a height. A seguir, abres o teu dataset usando a função open standard. Quando chamas o método read, não deixas os argumentos vazios. Argumentos vazios dizem ao Rasterio para carregar o ficheiro inteiro para a memória. Em vez disso, atribuis o teu novo objeto Window ao keyword argument window dentro do método read. O Rasterio traduz essa window em byte offsets no disco, vai buscar apenas esse chunk específico de 256 por 512, e devolve-o como um numpy array. O teu uso de memória mal se mexe. Se trabalhas com numpy o dia todo, passar offsets e dimensões pode parecer ligeiramente antinatural. Provavelmente estás habituado a definir slices com índices de start e stop para as tuas linhas e colunas. O Rasterio suporta este workflow com um método alternativo chamado Window dot from slices. Este método recebe dois inputs. O primeiro input é um tuple que define o start index e o stop index da linha. O segundo input é um tuple que define o start index e o stop index da coluna. Repara que a ordem aqui é linhas primeiro, e depois colunas. Isto espelha exatamente como fazes o slice de um numpy array bidimensional. Internamente, este método calcula os offsets e as dimensões por ti e devolve um objeto Window standard que podes passar ao método read tal como antes. Ao forçares windowed reads nos teus data pipelines, desacoplas permanentemente o tamanho físico dos teus ficheiros de input dos limites de hardware da tua máquina. É tudo por este episódio. Obrigado por ouvires, e continua a desenvolver!

Olá, daqui fala o Alex da DEV STORIES DOT EU. Rasterio: Numpy para Dados Geoespaciais, episódio 4 de 5. Acabaste de descarregar uma enorme imagem de satélite que cobre um estado inteiro, mas só te interessas pelos pixels dentro dos limites de uma quinta específica. Carregar a grid inteira desperdiça memória, e ler por pixel offsets é inútil quando o teu limite tem a forma de um polígono irregular. A solução é o Masking com vetores. Num episódio anterior, vimos o windowed reading. Essa abordagem exigia índices de pixel exatos, como pedir a row cinquenta e a column cem. O Masking é completamente diferente. Baseia-se em coordenadas geográficas do mundo real. Em vez de contar pixels, forneces uma vector geometry, como um polígono definido por longitudes e latitudes reais. O Rasterio descobre quais pixels calham dentro dessa forma e trata da tradução complexa entre o espaço contínuo de vector coordinates e a raster grid discreta. A ferramenta principal para este workflow é a função mask, localizada no módulo rasterio dot mask. Passas-lhe dois argumentos principais. Primeiro, um raster dataset aberto. Segundo, um iterable de formas geométricas. Estas formas devem ser formatadas como dictionaries que correspondam à especificação GeoJSON. Normalmente, lês um shapefile que contém o polígono da tua quinta, extrais o dictionary da geometry para essa feature específica, colocas tudo numa list padrão de Python, e passas à função mask. Aqui está o ponto chave. A função mask tem dois comportamentos distintos controlados por um único parâmetro boolean chamado crop. Precisas de perceber a diferença entre simplesmente fazer masking a um dataset e fazer-lhe efetivamente um crop. Se correres a função com o crop definido para false, o Rasterio olha para o polígono da tua quinta e define todos os pixels fora desse limite para um valor nodata. Ele lê o valor nodata padrão diretamente da metadata original do raster. O array de output que recebes tem exatamente as mesmas dimensões da tua enorme imagem original que cobre o estado. Os pixels dentro da quinta mantêm os seus valores reais, mas tudo o que está de fora é apagado. A extensão espacial geral do ficheiro não muda. Se definires o crop para true, o comportamento muda. O Rasterio continua a definir os pixels fora do teu polígono como nodata. No entanto, a seguir descarta as vastas áreas de espaço vazio. Ele calcula a bounding box retangular do polígono irregular da tua quinta e encolhe o array devolvido para caber exatamente dentro dessa bounding box específica. O teu gigante array estadual é fisicamente reduzido a uma grid minúscula e gerível, contendo apenas a quinta e uma pequena quantidade de padding de nodata à volta das suas extremidades irregulares. Como fazer crop altera fisicamente as dimensões da grid, também altera o ponto de partida geográfico. O canto superior esquerdo já não é o canto superior esquerdo do estado. É agora o canto superior esquerdo da bounding box da quinta. Portanto, a função mask devolve um tuple contendo dois itens. O primeiro item é o teu numpy array com o crop aplicado. O segundo item é um objeto de affine transform completamente novo. Este transform atualizado contém as novas coordenadas de origem. Se quiseres guardar este array com crop num novo ficheiro, tens de passar este transform atualizado para o teu write profile. Se reutilizares o transform original do estado por engano, a imagem da tua pequena quinta será esticada novamente por todo o estado. O verdadeiro valor da função mask é que ela trata de todo este recálculo espacial por ti automaticamente. Espero que tenha sido útil. Obrigado por ouvires, e aproveita o resto do teu dia.

Olá, daqui é o Alex da DEV STORIES DOT EU. Rasterio: Numpy para Dados Geoespaciais, episódio 5 de 5. Às vezes tens imagens de satélite belíssimas e de alta resolução, mas o teu modelo climático só aceita inputs mais grosseiros e de baixa resolução. Não podes simplesmente fazer slice do array, porque os pixels resultantes ainda precisam de representar exatamente a mesma área física, apenas numa grid de menor densidade. Para fazeres isso corretamente, precisas de Resampling e Writing. Isto não é sobre reprojection. O coordinate reference system permanece exatamente o mesmo. Estamos apenas a alterar a resolução. Digamos que precisas de fazer downsample de um raster para metade. No Rasterio, fazes o resampling no exato momento em que lês os dados para a memória. Quando chamas o método read no teu dataset, passas um parâmetro chamado out shape. Em vez de ler o array completo, o Rasterio lê-o e faz logo o scale para se ajustar ao shape que pediste. Para um raster de uma só band, o teu out shape seria um tuple a especificar uma band, depois metade da height original, e metade da width original. Também precisas de dizer ao Rasterio como calcular os valores destes novos pixels maiores. Fazes isso passando um parâmetro de resampling. O Rasterio fornece um enum Resampling para isto. Se estiveres a lidar com dados contínuos, como temperatura, podes passar Resampling dot bilinear, que calcula uma média ponderada dos pixels vizinhos. Se estiveres a lidar com dados categóricos, como land cover, usarias Resampling dot nearest para preservar os valores exatos das classes originais. Aqui está o ponto chave. Agora tens um novo array Numpy mais pequeno, mas que ainda representa exatamente a mesma área geográfica na Terra. Como o array tem menos pixels, cada pixel individual cobre mais espaço físico. Isto significa que o teu Affine transform original está agora completamente errado. Se escreveres o novo array para o disco usando o transform antigo, a tua imagem vai encolher no mapa para cobrir apenas um quarto da sua área original. Tens de fazer scale do transform. Fazes isto pegando no transform original do dataset e multiplicando-o por uma scaling matrix Affine. Encontras os scaling ratios dividindo a width original pela tua nova width, e a height original pela tua nova height. Como reduzimos as dimensões para metade, ambos estes ratios são exatamente dois. Quando multiplicas o transform original por estes scaling factors, obténs um novo transform onde as dimensões dos pixels são duplicadas, mas o ponto de origem no canto superior esquerdo permanece perfeitamente ancorado. Agora tens o teu array com resampling e o teu transform com scale. O passo final é escrever um GeoTIFF válido para o disco. Para fazeres isto, precisas de metadata. A abordagem mais segura é copiar a propriedade profile do teu dataset de origem. Isto dá-te um dictionary que contém o file driver, o coordinate reference system, o data type, e o valor nodata. Depois, atualizas este dictionary de profile com a tua nova height, a tua nova width, e o teu transform recém-calculado. Com o profile atualizado pronto, chamas rasterio dot open. Forneces o novo filename, defines o mode para write, e fazes unpack do teu dictionary de profile como keyword arguments. Isto cria um dataset vazio no disco, formatado de acordo com as tuas especificações exatas. Finalmente, chamas o método write neste novo objeto dataset e passas o teu array Numpy com downsample. O array é escrito para o ficheiro, e o teu novo GeoTIFF de menor resolução está pronto para o modelo climático. Lembra-te, modificar dados raster é sempre um contrato de duas partes: se alterares o shape do array Numpy, tens de alterar explicitamente o Affine transform para corresponder, ou o teu spatial footprint fica estragado. Isto leva-nos ao fim da nossa série. Encorajo-te a ler a documentação oficial do Rasterio para explorares outros algoritmos de resampling e tentares construir estas pipelines hands-on. Se quiseres sugerir tópicos para uma futura série, visita devstories dot eu. Gostaria de tirar um momento para te agradecer por ouvires — isso ajuda-nos muito. Tem um ótimo dia!