Grouping geommetries by colums

General Usage
, , ,
1

I need to combine geometries for a shapefile. Initially, there is shapefile_path_cities that will form the layer behind, after theres is shapefile_path_highways:

using GeoDataFrames
using DataFrames
using XLSX
using Plots
using Shapefile
using Random
using Statistics
using Base

###Shapefile municipalities

shapefile_path = "C:\\Users\\daves\\OneDrive\\Pessoal\\Acadêmico\\Mestrado\\Dissertação - Execução\\Análises\\MT_Municipios_2022\\MT_Municipios_2022.shp"
municipalities = Shapefile.Table(shapefile_path) |> DataFrame

# Highway shapefile
shapefile_path_highways = "C:\\Users\\daves\\OneDrive\\Pessoal\\Acadêmico\\Mestrado\\Dissertação - Execução\\Análises\\Rodovias\\SNV_202410A.shp"
highways = Shapefile.Table(shapefile_path_rodovias) |> DataFrame

These data on shapefile_path_highway, which be puted in the dataframe highways will be filtered and grouped:

# Filtering sg_uf equals to MT
filtered_highwayss_MT = filter(row -> row.sg_uf == "MT", rodovias)

# Grouping by vl_br row
grouped_highways = groupby(filtered_highways_MT, :vl_br)

Untill now everthing done succesfully, however the combine function did not work well:

#Combine geometries
function combine_geometries(geometries)
    combined_geom = ArchGDAL.readgeometry(geometries[1].geometry)
    for geom in geometries[2:end]
        combined_geom = ArchGDAL.union!(combined_geom, ArchGDAL.readgeometry(geom.geometry))
    end
    return combined_geom
end

# Applying function
combined_geometries = combine(grouped_rodovias, :geometry => combine_geometries => :geometry)

#Showing results
println(combined_geometries)

If you want to take shapefiles the first is in and https://geoftp.ibge.gov.br/organizacao_do_territorio/malhas_territoriais/malhas_municipais/municipio_2022/UFs/MT/MT_Municipios_2022.zip the second is in https://servicos.dnit.gov.br/dnitcloud/index.php/s/oTpPRmYs5AAdiNr/download. The first is provided by IBGE and the second by DNIT, both agencies of Brazilian government.

There is an argument to combine geometries that can function instead of readgeometry and union!?

2

Shapefile.jl geometries are GeoInterface compatible, so you should use GeoInterface.convert(ArchGDAL, geometry) to convert them. You will likely get an error message trying to use ArchGDAL internal functions on Shapefile.jl geometries (also I can’t find any function named readgeometry in ArchGDAL? union! is also not defined for ArchGDAL geometries…)

Try something like this:

function combine_geometries(geometries)

    return mapfoldl(
        x -> GeoInterface.convert(ArchGDAL, x), # map - apply to each element
        ArchGDAL.union, # reduce 2 elements to 1 element
        geometries # the vector to perform the operation on
    )
end
1 Like
3

Thanks, this code worked well!