Módulo 02

Objetivo

A partir de um arquivo vetorial (shape), realizar a reprojeção para outro CRS (DATUM).

Leitura de arquivo (shape) e reprojeção

Importação dos pacotes

Inicialmente importar os pacotes

from osgeo import ogr, osr, gdal

Leitura do arquivo vetorial

Para abrir o arquivo vetorial, define-se o caminho e o parâmetro 0 ou 1:

# reading shape
shapePath = './shape/mt_municipios_2020.shp'
ds = ogr.Open(shapePath, 1)  # 0 means read-only. 1 means writeable

# Check to see if shapefile was found.
if ds is None:
    print("------------------------------------------------------- \n")
    print('Could not open %s' % (shapePath))
    print('Check shapefile')
else:
    print("------------------------------------------------------- \n")
    print('%s' ' has been opened' % (shapePath))
    layer = ds.GetLayer()

Nesse bloco foi especificado o caminho que o arquivo encontra-se, e o método Open da classe OGR é executado com o argumento 1 para possibilitar a edição do arquivo. Na sequência, criou-se uma verificação da leitura do arquivo. Se o arquivo foi aberto com sucesso, o objeto layer recebe a classe Layer via ds.GetLayer(). Verifique aqui todos os métodos disponíveis na classe Layer aqui

Reprojeção do arquivo

Como já comentado, é possível GDAL abrir mais de 200 tipos de arquivos vetoriais e raster, para saber quais são os drivers disponíveis é possível usar o seguinte código:

# ------------------------------------------------------------
# Drivers in gdal
# ------------------------------------------------------------
print(gdal.GetDriverCount())
lsDrivers = [(i, gdal.GetDriver(i).GetDescription())
             for i in range(gdal.GetDriverCount())]
print(lsDrivers)

> [(0, 'VRT'), (1, 'DERIVED'), (2, 'GTiff'), (3, 'NITF'), (4, 'RPFTOC'), (5, 'ECRGTOC'), (6, 'HFA'), (7, 'SAR_CEOS'), (8, 'CEOS'), (9, 'JAXAPALSAR'), (10, 'GFF'), (11, 'ELAS'), (12, 'AIG'), (13, 'AAIGrid'), (14, 'GRASSASCIIGrid'), (15, 'SDTS'), (16, 'DTED'), (17, 'PNG'), (18, 'JPEG'), (19, 'MEM'), (20, 'JDEM'), (21, 'GIF'), (22, 'BIGGIF'), (23, 'ESAT'), (24, 'BSB'), (25, 'XPM'), (26, 'BMP'), (27, 'DIMAP'), (28, 'AirSAR'), (29, 'RS2'), (30, 'SAFE'), (31, 'PCIDSK'), (32, 'PCRaster'), (33, 'ILWIS'), (34, 'SGI'), (35, 'SRTMHGT'), (36, 'Leveller'), (37, 'Terragen'), (38, 'GMT'), (39, 'netCDF'), (40, 'HDF4'), (41, 'HDF4Image'), (42, 'ISIS3'), (43, 'ISIS2'), (44, 'PDS'), (45, 'PDS4'), (46, 'VICAR'),... (110, 'BAG'), (111, 'HDF5'), (112, 'HDF5Image'), (113, 'NWT_GRD'), (114, 'NWT_GRC'), (115, 'ADRG'), (116, 'SRP'), (117, 'BLX'), (118, 'EPSILON'), (119, 'PostGISRaster'), (120, 'SAGA'), (121, 'XYZ'), (122, 'HF2'), (123, 'JPEGLS'), (124, 'OZI'), (125, 'CTG'), (126, 'E00GRID'), (127, 'ZMap'), (128, 'NGSGEOID'), (129, 'IRIS'), (130, 'PRF'), (131, 'RDA'), (132, 'EEDAI'), (133, 'EEDA'), (134, 'SIGDEM'), (135, 'IGNFHeightASCIIGrid'), (136, 'GNMFile'), (137, 'GNMDatabase'), (138, 'ESRI Shapefile'), ..., (221, 'EHdr'), (222, 'ISCE'), (223, 'HTTP')]

Uma vez que já conhecemos os possíveis tipos de arquivo para fazer a reprojeção, a primeira etapa é saber o CRS do arquivo:

# get CRS
>>> print(layer.GetSpatialRef())
GEOGCS["SIRGAS 2000",
    DATUM["Sistema_de_Referencia_Geocentrico_para_las_AmericaS_2000",
        SPHEROID["GRS 1980",6378137,298.257222101,
            AUTHORITY["EPSG","7019"]],
        TOWGS84[0,0,0,0,0,0,0],
        AUTHORITY["EPSG","6674"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4674"]]

Conforme o código acima, verifica-se que o CRS do arquivo em questão é o EPSG:4674 (“SIRGAS 2000”). Neste exercício faremos a reprojeção para o CRS EPSG:5880 (“SIRGAS 2000 / Brazil Polyconic”).

# create a sourceCrs and a targetCrs
>>> sourceCrs = osr.SpatialReference()
>>> sourceCrs.ImportFromEPSG(4674)
0
>>> targetCrs = osr.SpatialReference()
>>> targetCrs.ImportFromEPSG(5880)
0

Agora que o CRS do arquivo de entrada é conhecido, e os objetos de referência dos CRS foram criados, é possível preparar o arquivo vetorial (shape):

# assign driver
driver = ogr.GetDriverByName('ESRI Shapefile')

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(sourceCrs, targetCrs)

# create the output layer as ogr.wkbMultiPolygon
outputShapefile = './shape/mt_municipios_2020_crs5880.shp'
if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)
outDs = driver.CreateDataSource(outputShapefile)
outLayer = outDs.CreateLayer(
    "mt_municipios_2020_crs5880", targetCrs, geom_type=ogr.wkbMultiPolygon)

# add fields and defn from layer (original)
inLayerDefn = layer.GetLayerDefn()
for idx in range(inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(idx)
    outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
inFeature = layer.GetNextFeature()
while inFeature:
    # get the input geometry
    geom = inFeature.GetGeometryRef()
    # reproject the geometry
    geom.Transform(coordTrans)
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(
            i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    outLayer.CreateFeature(outFeature)
    # dereference the features and get the next input feature
    outFeature = None
    inFeature = layer.GetNextFeature()

# Save and close the shapefiles
inDataSet = outDs = outFeature = inFeature = None

The power of GDAL and OGR

Leitura de arquivo (shape) e reprojeção (CRS)