import geopandas as gpd
# open the datasets
brazil = gpd.read_file('brazil')
stations = gpd.read_file('ozone_stations.csv')
# create the geometry for the stations
stations['geometry'] = gpd.points_from_xy(x=stations['X'], y=stations['Y'], crs='epsg:4326')
# filter just the stations in Brazil
br_stations = stations.query("country == 'Brazil'")
# create the coverage area for each station
br_stations.geometry = br_stations.buffer(5)
# project both datasets to equal area CRS
brazil = brazil.to_crs('ESRI:54034')
br_stations = br_stations.to_crs('ESRI:54034')
# perform an overlay difference
uncovered = brazil.overlay(br_stations, how='difference')
# calculate the not covered area
ncarea = format(round(uncovered.area.sum()/1e6), ',')
# plot the results
ax = uncovered.plot(facecolor='green', edgecolor='black')
ax.set_ylabel('Latitude (deg)')
ax.set_xlabel('Longitude (deg)')
ax.set_title(f'Uncovered area is {ncarea} km²')
ax.figure.savefig('output/uncovered.png', dpi=300)