Recorte e plot de shapefile

Um shapefile é um formato de armazenamento de dados de vetor da Esri para armazenar a posição, a forma e os atributos de feições geográficas. Uma das formas de visualizar um shapefile online é através do site mapshaper. Outra forma é usar um script simples em python, apresentado nesse post. Como o arquivo original é muito grande, primeiro ele será recortado usando um comando no terminal Linux.

Shapefile recortado e plotado seguindo os procedimentos do post
Shapefile recortado e plotado seguindo os procedimentos do post

O arquivo original compreende todas as massas d’água do país, podendo ser baixado no formato SHP no Catálogo de Metadados da ANA no link. Como o arquivo é muito grande e só traablharemos com uma região pequena dele, o recorte pode ser feito usando o programa GDAL através do seguinte comando na pasta do arquivo baixado:

ogr2ogr -clipsrc -63 -6 -57 0 Amazon_Rivers.shp geoft_bho_massa_dagua_v2019.shp

O parâmetro “clipsrc” executa o recorte entre as longitudes (x) e latitudes (y) informadas na seguinte ordem: [xmin ymin xmax ymax]. Diferentemente da maioria dos comandos, o primeiro arquivo informado é o de saída e o segundo, o de entrada – ou seja, o arquivo gerado tem o nome “Amazon_Rivers.shp”.

O script a seguir usa os pacotes cartopy e matplotlib para gerar uma imagem do shapefile. Um “mapa base” é feito na mesma área recortada com os contornos continentais (particularmente nessa imagem não são necessários). Na sequência, é colocada cor nas superfícies de terra e então é lido e plotado o shapefile. A imagem finalmente é ajustada e impressa para um arquivo do tipo PNG.

#!/usr/bin/env python3.7.11
# -*- Coding: UTF-8 -*-

import cartopy.crs as ccrs
import cartopy.feature as cf
import cartopy.io.shapereader as shpreader
import cartopy.feature as cfeature
from matplotlib import pyplot as plt

# Plot base map
ax = plt.axes(projection=ccrs.PlateCarree())
ax.set_extent([-61,-59,-4,-2])
# Show continents lines
ax.add_feature(cf.COASTLINE)

# Display land shape
land_feature = cfeature.NaturalEarthFeature('physical', 'land',\
 '110m', edgecolor='face', facecolor=cfeature.COLORS['land'])
ax.add_feature(land_feature)

# Read shape file
reader = shpreader.Reader('helpers/Amazon_Rivers.shp')
# Display shape
shape_feature = cfeature.ShapelyFeature(reader.geometries(),\
 ccrs.PlateCarree(), facecolor=cfeature.COLORS['water'], edgecolor='black', lw=0.4)
ax.add_feature(shape_feature)

# Tight layout to reduce blank/empty space around map
plt.tight_layout()
# Save figure as SVG
plt.savefig('map.png')

Fontes

Compartilhe :)

Leave a Reply

O seu endereço de e-mail não será publicado.

Esse site utiliza o Akismet para reduzir spam. Aprenda como seus dados de comentários são processados.