# 1. Project Setup
Initially, we want to load libraries. Be certain that all of them are properly installed. Also, I’m using Python 3.12.3.
## 1.1 Load libraries# If it's essential install any library, run below:
# pip install library1 library2 library3 ...
# Basic Libraries
import os # For file operations
import gc # For garbage collection, it avoids RAM memory issues
import numpy as np # For coping with matrices
import pandas as pd # For coping with dataframes
import janitor # For data cleansing (mainly column names)
import numexpr # For fast pd.query() manipulation
import inflection # For string manipulation
import unidecode # For string manipulation
# Geospatial Libraries
import geopandas as gpd # For coping with shapefiles
import pyogrio # For fast .gpkg file manipulation
import ee # For Google Earth Engine API
import contextily as ctx # For basemaps
import folium # For interactive maps
# Shapely Objects and Geometry Manipulation
from shapely.geometry import mapping, Polygon, Point, MultiPolygon, LineString # For geometry manipulation
# Raster Data Manipulation and Visualization
import rasterio # For raster manipulation
from rasterio.mask import mask # For raster data manipulation
from rasterio.plot import show # For raster data visualization
# Plotting and Visualization
import matplotlib.pyplot as plt # For plotting and data visualization
from matplotlib.colours import ListedColormap, Normalize # For color manipulation
import matplotlib.colours as colours # For color manipulation
import matplotlib.patches as mpatches # For creating patch objects
import matplotlib.cm as cm # For colormaps
# Data Storage and Manipulation
import pyarrow # For efficient data storage and manipulation
# Video Making
from moviepy.editor import ImageSequenceClip # For creating videos (section 4.7 only) - check this if you may have issues: https://github.com/kkroening/ffmpeg-python
Then, be certain that you may have a folder for this project. All resources and outputs might be saved there. This folder could be situated in your local drive, a cloud-based storage solution, or in a selected folder on Google Drive where you’ll save the rasters retrieved using the GEE API.
When running your code, be certain that to vary the address below to your project path. Windows users, all the time remember to make use of as an alternative of /.
# 1.2 Set working directory
project_path = 'path_to_your_project_folder' # Where you'll save all outcomes and resources have to be in
os.chdir(project_path) # All resources on the project are relative to this path# 1.3 Further settings
pd.set_option('compute.use_numexpr', True) # Use numexpr for fast pd.query() manipulation
Lastly, this function is beneficial for plotting geometries over OpenStreetMap (OSM). It is especially helpful when working with unknown shapefiles to make sure accuracy and avoid mistakes.
## 1.4 Set function to plot geometries over an OSM
def plot_geometries_on_osm(geometries, zoom_start=10):# Calculate the centroid of all geometries to center the map
centroids = [geometry.centroid for geometry in geometries]
avg_x = sum(centroid.x for centroid in centroids) / len(centroids)
avg_y = sum(centroid.y for centroid in centroids) / len(centroids)
# Create a folium map centered around the common centroid
map = folium.Map(location=[avg_y, avg_x], zoom_start=zoom_start)
# Add each geometry to the map
for geometry in geometries:
geojson = mapping(geometry) # Convert the geometry to GeoJSON
folium.GeoJson(geojson).add_to(map)
return map
# 2. Single Example: Acrelândia (AC) in 2022
For example to create intuition of the method, we’ll save, clean, and plot land use in Acrelândia (AC) in 2022. It’s a city in the midst of the AMACRO region (the three-state border of Amazonas, Acre, and Rondônia), where the usually untouched forest is being rapidly destroyed.
On this section, I’ll explain step-by-step of the script, after which standardize the method to run it for multiple places over multiple years. Since saving large rasters using the API generally is a slow process, I like to recommend using it provided that it’s essential take care of just a few or small areas for just a few years. Large areas may take hours to avoid wasting on Google Drive, so I like to recommend downloading the heavy LULC files for the entire country after which cleansing them, as we’ll do in a future post.
To run the code, first download and save IBGE’s¹ Brazilian cities shapefiles (select Brazil > Municipalities). Remember, you need to use any shapefile in Brazil to perform this algorithm.
## 2.1 Get the geometry of the realm of interest (Acrelândia, AC)
brazilian_municipalities = gpd.read_file('municipios/file.shp', engine='pyogrio', use_arrow=True) # Read the shapefile - you need to use every other shapefile here. Shapefiles have to be in your project folder, as set in 1.2
brazilian_municipalities = brazilian_municipalities.clean_names() # Clean the column names (remove special characters, spaces, etc.)
brazilian_municipalities.crs = 'EPSG:4326' # Set the CRS to WGS84 (MapBiomas uses this CRS)
brazilian_municipalities
## 2.2 Get geometry for Acrelândia, AC
city = brazilian_municipalities.query('nm_mun == "Acrelândia"') # Filter the geometry for Acrelândia, AC (could be every other city or set of cities)
city_geom = city.geometry.iloc[0] # Get the geometry of Acrelândia, AC
city_geom # See the geometry shape
Once we’ve got the shapefile we wish to check properly saved, we’ll create a bounding box around it to crop the MapBiomas full raster. Then, we’ll reserve it the GEE Python API.
## 2.3 Set the bounding box (bbox) for the realm of interest
bbox = city_geom.bounds # Get the bounding box of the geometry
bbox = Polygon([(bbox[0], bbox[1]), (bbox[0], bbox[3]), (bbox[2], bbox[3]), (bbox[2], bbox[1])]) # Convert the bounding box to a Polygonbbox_xmin = bbox.bounds[0] # Get the minimum x coordinate of the bounding box
bbox_ymin = bbox.bounds[1] # Get the minimum y coordinate of the bounding box
bbox_xmax = bbox.bounds[2] # Get the utmost x coordinate of the bounding box
bbox_ymax = bbox.bounds[3] # Get the utmost y coordinate of the bounding box
bbox # See bbox around Acrelândia shape
# Plot the bounding box and the geometry of Acrelandia over an OSM map
plot_geometries_on_osm([bbox, city_geom], zoom_start=10)
Now, we’ll access the MapBiomas Google Earth Engine API. First, we want to create a cloud project on GEE using a Google Account. Be certain that you may have enough space in your Google Drive account.
Then, we want to authenticate the GEE Python API (just once). In case you are a VSCode user, notice that the token insertion box appears in the highest right corner of the IDE.
All images from the MapBiomas LULC Collection can be found in the identical asset. Notice that you could barely modify this script to work with other assets within the GEE catalog and other MapBiomas collections.
## 2.4 Acess MapBiomas Collection 8.0 using GEE API
# import ee - already imported at 1.1ee.Authenticate() # Just for the primary time
ee.Initialize() # Run it each time you begin a brand new session
# Define the MapBiomas Collection 8.0 asset ID - retrieved from https://brasil.mapbiomas.org/en/colecoes-mapbiomas/
mapbiomas_asset = 'projects/mapbiomas-workspace/public/collection8/mapbiomas_collection80_integration_v1'
asset_properties = ee.data.getAsset(mapbiomas_asset) # Check the asset's properties
asset_properties # See properties
Here, each band represents the LULC data for a given yr. Be certain that that the code below is correctly written. This selects the image for the specified yr after which crops the raw raster for a bounding box across the region of interest (ROI) — Acrelândia, AC.
## 2.5 Filter the gathering for 2022 and crop the gathering to a bbox around Acrelândia, AC
yr = 2022
band_id = f'classification_{yr}' # bands (or yearly rasters) are named as classification_1985, classification_1986, ..., classification_2022mapbiomas_image = ee.Image(mapbiomas_asset) # Get the photographs of MapBiomas Collection 8.0
mapbiomas2022 = mapbiomas_image.select(band_id) # Select the image for 2022
roi = ee.Geometry.Rectangle([bbox_xmin, bbox_ymin, bbox_xmax, bbox_ymax]) # Set the Region of Interest (ROI) to the bbox around Acrelândia, AC - set in 2.3
image_roi = mapbiomas2022.clip(roi) # Crop the image to the ROI
Now, we save the cropped raster on Google Drive (in my case, into the ‘tutorial_mapbiomas_gee’ folder). Be certain that you may have created the destination folder in your drive before running.
I attempted to reserve it locally, nevertheless it looks like it’s essential save GEE rasters at Google Drive (let me know should you know tips on how to do it locally). That is essentially the most time-consuming a part of the code. For big ROIs, this might take hours. Check your GEE task manager to see if the rasters were properly loaded to the destination folder.
## 2.6 Export it to your Google Drive (ensure you may have space there and that it is correctly arrange)# Obs 1: Recall it's essential authenticate the session, because it was done on 2.4
# Obs 2: Ensure you may have enough space on Google Drive. Free version only gives 15 Gb of storage.
export_task = ee.batch.Export.image.toDrive(
image=image_roi, # Image to export to Google Drive as a GeoTIFF
description='clipped_mapbiomas_collection8_acrelandia_ac_2022', # Task description
folder='tutorial_mapbiomas_gee', # Change this to the folder in your Google Drive where you wish to save the file
fileNamePrefix='acrelandia_ac_2022', # File name (change it if you wish to)
region=roi.getInfo()['coordinates'], # Region to export the image
scale=30,
fileFormat='GeoTIFF'
)
# Start the export task
export_task.start()
# 3. Plot the Map
Now we’ve got a raster with LULC data for a bounding box around Acrelândia in 2022. That is saved on the address below (at Google Drive). First, let’s see the way it looks.
## 3.1 Plot the orginal raster over a OSM
file_path = 'path_of_exported_file_at_google_drive.tif' # Change this to the trail of the exported file# Plot data
with rasterio.open(file_path) as src:
data = src.read(1)
print(src.meta)
print(src.crs)
show(data)
In MapBiomas LULC Collection 8, each pixel represents a selected land use type in accordance with this list. As an illustration, ‘3’ means ‘Natural Forest’, ‘15’ means ‘Pasture’, and ‘0’ means ‘No data’ (pixels within the raster not inside the Brazilian borders).
We are going to explore the info we’ve got before plotting it.
## 3.2 See unique values
unique_values = np.unique(data)
unique_values # Returns unique pixels values within the raster# 0 = no data, parts of the image outside Brazil
## 3.3 See the frequency of every class (except 0 - no data)
unique_values, counts = np.unique(data[data != 0], return_counts=True) # Get the unique values and their counts (except zero)
pixel_counts = pd.DataFrame({'value': unique_values, 'count': counts})
pixel_counts['share'] = (pixel_counts['count'] / pixel_counts['count'].sum())*100
pixel_counts
At the top of the day, we’re working with a big matrix wherein each element represents how each tiny 30m x 30m piece of land is used.
## 3.4 See the actual raster (a matrix wherein each element represents a pixel value - land use code on this case)
data
Now, we want to prepare our raster data. As an alternative of categorizing each pixel by exact land use, we’ll categorize them more broadly. We are going to divide pixels into natural forest, natural non-forest vegetation, water, pasture, agriculture, and other uses. Specifically, we’re concerned about tracking the conversion of natural forest to pasture. To realize this, we’ll reassign pixel values based on the mapbiomas_categories dict below, which follows with MapBiomas’ land use and land cover (LULC) categorization.
The code below crops the raster to Acrelândia’s limits and reassigns pixels in accordance with the mapbiomas_categories dict. Then, it saves it as a brand new raster at ‘reassigned_raster_path’. Note that the old raster was saved on Google Drive (after being downloaded using GEE’s API), while the brand new one might be saved within the project folder (in my case, a OneDrive folder on my PC, as set in section 1.2). From here onwards, we’ll use only the reassigned raster to plot the info.
That is the important a part of the script. If you may have doubts about what is going on here (cropping for Acrelândia after which reassigning pixels to broader categories), I like to recommend running it and printing results for each step.
mapbiomas_categories = {
# Forest (= 3)
1:3, 3:3, 4:3, 5:3, 6:3, 49:3, # That's, values 1, 3, 4, 5, 6, and 49 might be reassigned to three (Forest)
# Other Non-Forest Natural Vegetation (= 10)
10:10, 11:10, 12:10, 32:10, 29:10, 50:10, 13:10, # That's, values 10, 11, 12, 32, 29, 50, and 13 might be reassigned to 10 (Other Non-Forest Natural Vegetation)
# Pasture (= 15)
15:15,
# Agriculture (= 18)
18:18, 19:18, 39:18, 20:18, 40:18, 62:18, 41:18, 36:18, 46:18, 47:18, 35:18, 48:18, 21:18, 14:18, 9:18, # That's, values 18, 19, 39, 20, 40, 62, 41, 36, 46, 47, 35, 48, 21, 14, and 9 might be reassigned to 18 (Agriculture)
# Water ( = 26)
26:26, 33:26, 31:26, # That's, values 26, 33, and 31 might be reassigned to 26 (Water)
# Other (= 22)
22:22, 23:22, 24:22, 30:22, 25:22, 27:22, # That's, values 22, 23, 24, 30, 25, and 27 might be reassigned to 22 (Other)
# No data (= 255)
0:255 # That's, values 0 might be reassigned to 255 (No data)
}
## 3.5 Reassing pixels values to the MapBiomas custom general categories and crop it to Acrelandia, AC limits
original_raster_path = 'path_to_your_google_drive/tutorial_mapbiomas_gee/acrelandia_ac_2022.tif'
reassigned_raster_path = 'path_to_reassigned_raster_at_project_folder' # Somewhere within the project folder set at 1.2with rasterio.open(original_raster_path) as src:
raster_array = src.read(1)
out_meta = src.meta.copy() # Get metadata from the unique raster
# 3.5.1. Crop (or mask) the raster to the geometry of city_geom (on this case, Acrelandia, AC) and thus remove pixels outside the town limits (might be assigned to no data = 255)
out_image, out_transform = rasterio.mask.mask(src, [city_geom], crop=True)
out_meta.update({
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform
}) # Update metadata to the brand new raster
raster_array = out_image[0] # Get the masked raster
modified_raster = np.zeros_like(raster_array) # Base raster filled with zeros to be modified
# 3.5.2. Reassign each pixel based on the mapbiomas_categories dictionary
for original_value, new_value in mapbiomas_categories.items():
mask = (raster_array == original_value) # Create a boolean mask for the unique value (True = Replace, False = Don't replace)
modified_raster[mask] = new_value # Replace the unique values with the brand new values, when needed (that's, when the mask is True)
out_meta = src.meta.copy() # Get metadata from the unique raster
out_meta.update(dtype=rasterio.uint8, count=1) # Update metadata to the brand new raster
with rasterio.open(reassigned_raster_path, 'w', **out_meta) as dst: # Write the modified raster to a brand new file on the reassigned_raster_path
dst.write(modified_raster.astype(rasterio.uint8), 1)
## 3.6 See the frequency of pixels within the reassigned raster
with rasterio.open(reassigned_raster_path) as src:
raster_data = src.read(1)
unique_values = np.unique(raster_data)
total_non_zero = np.sum(raster_data != 255) # Count the entire variety of non-zero pixelsfor value in unique_values:
if value != 255: # Exclude no data (= 255)
count = np.sum(raster_data == value) # Count the variety of pixels with the worth
share = count / total_non_zero # Calculate the share of the worth
share = share.round(3) # Round to three decimal places
print(f"Value: {value}, Count: {count}, Share: {share}")
Now we plot the info with generic colours. We are going to enhance the map later, but that is just a primary (or second?) look. Notice that I specifically set 255 (= no data, pixels outside Acrelândia) to be white for higher visualization.
## 3.7 Plot the reassigned raster with generic colours
with rasterio.open(reassigned_raster_path) as src:
data = src.read(1) # Read the raster data
unique_values = np.unique(data) # Get the unique values within the rasterplt.figure(figsize=(10, 8)) # Set the figure size
cmap = plt.cm.viridis # Using Viridis colormap
norm = Normalize(vmin=data.min(), vmax=26) # Normalize the info to the range of the colormap (max = 26, water)
masked_data = np.ma.masked_where(data == 255, data) # Mask no data values (255)
plt.imshow(masked_data, cmap=cmap, norm=norm) # Plot the info with the colormap
plt.colorbar(label='Value') # Add a colorbar with the values
plt.show()
Now it’s time to create an attractive map. I actually have chosen Matplotlib because I need static maps. In case you prefer interactive choropleths, you need to use Plotly.
For further details on choropleths with Matplotlib, check its documentation, GeoPandas guide, and the good Yan Holtz’s Python Graph Gallery — where I get much of the inspiration and tools for DataViz in Python. Also, for beautiful color palettes, coolors.co is a wonderful resource.
Be certain that you may have all data visualization libraries properly loaded to run the code below. I also tried to vary the order of patches, but I didn’t know tips on how to. Let me know should you discover tips on how to do it.
## 3.8 Plot the reassigned raster with custom colours# Define the colours for every class - notice it's essential follow the identical order because the values and have to be numerically increasing or decreasing (still need to search out out tips on how to solve it)
values = [3, 10, 15, 18, 22, 26, 255] # Values to be coloured
colors_list = ['#6a994e', '#a7c957', '#c32f27', '#dda15e', '#6c757d', '#0077b6','#FFFFFF'] # HEX codes of the colours used
labels = ['Natural Forest', 'Other Natural Vegetation', 'Pasture', 'Agriculture', 'Others', 'Water', 'No data'] # Labels displayed on the legend
cmap = colours.ListedColormap(colors_list) # Create a colormap (cmap) with the colours
bounds = values + [256] # Add a worth to the top of the list to incorporate the last color
norm = colours.BoundaryNorm(bounds, cmap.N) # Normalize the colormap to the values
img = plt.imshow(raster_data, interpolation='nearest', cmap=cmap, norm=norm) # Plot the info with the colormap
legend_patches = [mpatches.Patch(color=colors_list[i], label=labels[i]) for i in range(len(values)-1)] # Create the legend patches withou the last one (255 = no data)
# Create the legend
plt.legend(handles = legend_patches, # Add the legend patches
bbox_to_anchor = (0.5, -0.02), # Place the legend below the plot
loc = 'upper center', # Place the legend within the upper center
ncol = 3, # Variety of columns
fontsize = 9, # Font size
handlelength=1,# Length of the legend handles
frameon=False) # Remove the frame across the legend
plt.axis('off') # Remove the axis
plt.title('Land Use in Acrelândia, AC (2022)', fontsize=20) # Add title
plt.savefig('figures/acrelandia_ac_2022.pdf', format='pdf', dpi=1800) # Reserve it as a PDF on the figures folder
plt.show()
4. Standardized Functions
Now that we’ve got built intuition on tips on how to download, save, clean, and plot MapBiomas LULC rasters. It’s time to generalize the method.
On this section, we’ll define functions to automate these steps for any shape and any yr. Then, we’ll execute these functions in a loop to research a selected city — Porto Acre, AC — from 1985 to 2022. Finally, we’ll make a video illustrating the LULC evolution in that area over the desired period.
First, save a bounding box (bbox) across the Region of Interest (ROI). You simply must input the specified geometry and specify the yr. This function will save the bbox raster across the ROI to your Google Drive.
## 4.1 For a generic geometry in any yr, save a bbox across the geometry to Google Drivedef get_mapbiomas_lulc_raster(geom, geom_name, yr, folder_in_google_drive):
ee.Authenticate() # Just for the primary time
ee.Initialize() # Run it each time you begin a brand new session
my_geom = geom
bbox = my_geom.bounds # Get the bounding box of the geometry
bbox = Polygon([(bbox[0], bbox[1]), (bbox[0], bbox[3]), (bbox[2], bbox[3]), (bbox[2], bbox[1])]) # Convert the bounding box to a Polygon
bbox_xmin = bbox.bounds[0] # Get the minimum x coordinate of the bounding box
bbox_ymin = bbox.bounds[1] # Get the minimum y coordinate of the bounding box
bbox_xmax = bbox.bounds[2] # Get the utmost x coordinate of the bounding box
bbox_ymax = bbox.bounds[3] # Get the utmost y coordinate of the bounding box
mapbiomas_asset = 'projects/mapbiomas-workspace/public/collection8/mapbiomas_collection80_integration_v1'
band_id = f'classification_{yr}'
mapbiomas_image = ee.Image(mapbiomas_asset) # Get the photographs of MapBiomas Collection 8.0
mapbiomas_data = mapbiomas_image.select(band_id) # Select the image for 2022
roi = ee.Geometry.Rectangle([bbox_xmin, bbox_ymin, bbox_xmax, bbox_ymax]) # Set the Region of Interest (ROI) to the bbox around the specified geometry
image_roi = mapbiomas_data.clip(roi) # Crop the image to the ROI
export_task = ee.batch.Export.image.toDrive(
image=image_roi, # Image to export to Google Drive as a GeoTIFF
description=f"save_bbox_around_{geom_name}_in_{yr}", # Task description
folder=folder_in_google_drive, # Change this to the folder in your Google Drive where you wish to save the file
fileNamePrefix=f"{geom_name}_{yr}", # File name
region=roi.getInfo()['coordinates'], # Region to export the image
scale=30,
fileFormat='GeoTIFF'
)
export_task.start() # Notice that uploading those rasters to Google Drive may take some time, specially for giant areas
# Test it using Rio de Janeiro in 2022
folder_in_google_drive = 'tutorial_mapbiomas_gee'
rio_de_janeiro = brazilian_municipalities.query('nm_mun == "Rio de Janeiro"')
rio_de_janeiro.crs = 'EPSG:4326' # Set the CRS to WGS84 (this project default one, change if needed)
rio_de_janeiro_geom = rio_de_janeiro.geometry.iloc[0] # Get the geometry of Rio de Janeiro, RJteste1 = get_mapbiomas_lulc_raster(rio_de_janeiro_geom, 'rio_de_janeiro', 2022, folder_in_google_drive)
Second, crop the raster to incorporate only the pixels inside the geometry and reserve it as a brand new raster.
I selected to reserve it on Google Drive, but you’ll be able to change `reassigned_raster_path` to reserve it anywhere else. In case you change it, be certain that to update the remainder of the code accordingly.
Also, you’ll be able to reassign pixels as needed by modifying the mapbiomas_categories dict. The left digit represents the unique pixel values, and the best one represents the reassigned (recent) values.
## 4.2 Crop the raster for the specified geometry
def crop_mapbiomas_lulc_raster(geom, geom_name, yr, folder_in_google_drive):
original_raster_path = f'path_to_your_google_drive/{folder_in_google_drive}/{geom_name}_{yr}.tif'
reassigned_raster_path = f'path_to_your_google_drive/{folder_in_google_drive}/cropped_{geom_name}_{yr}.tif'my_geom = geom
mapbiomas_categories = {
# Forest (= 3)
1:3, 3:3, 4:3, 5:3, 6:3, 49:3,
# Other Non-Forest Natural Vegetation (= 10)
10:10, 11:10, 12:10, 32:10, 29:10, 50:10, 13:10,
# Pasture (= 15)
15:15,
# Agriculture (= 18)
18:18, 19:18, 39:18, 20:18, 40:18, 62:18, 41:18, 36:18, 46:18, 47:18, 35:18, 48:18, 21:18, 14:18, 9:18,
# Water ( = 26)
26:26, 33:26, 31:26,
# Other (= 22)
22:22, 23:22, 24:22, 30:22, 25:22, 27:22,
# No data (= 255)
0:255
} # You may change this to whatever categorization you wish, but just remember to adapt the colours and labels within the plot
with rasterio.open(original_raster_path) as src:
raster_array = src.read(1)
out_meta = src.meta.copy() # Get metadata from the unique raster
# Crop the raster to the geometry of my_geom and thus remove pixels outside the town limits (might be assigned to no data = 0)
out_image, out_transform = rasterio.mask.mask(src, [my_geom], crop=True)
out_meta.update({
"height": out_image.shape[1],
"width": out_image.shape[2],
"transform": out_transform
}) # Update metadata to the brand new raster
raster_array = out_image[0] # Get the masked raster
modified_raster = np.zeros_like(raster_array) # Base raster filled with zeros to be modified
# Reassign each pixel based on the mapbiomas_categories dictionary
for original_value, new_value in mapbiomas_categories.items():
mask = (raster_array == original_value) # Create a boolean mask for the unique value (True = Replace, False = Don't replace)
modified_raster[mask] = new_value # Replace the unique values with the brand new values, when needed (that's, when the mask is True)
out_meta = src.meta.copy() # Get metadata from the unique raster
out_meta.update(dtype=rasterio.uint8, count=1) # Update metadata to the brand new raster
with rasterio.open(reassigned_raster_path, 'w', **out_meta) as dst: # Write the modified raster to a brand new file on the reassigned_raster_path
dst.write(modified_raster.astype(rasterio.uint8), 1)
teste2 = crop_mapbiomas_lulc_raster(rio_de_janeiro_geom, 'rio_de_janeiro', 2022, folder_in_google_drive)
Now we see the frequency of every pixel within the cropped reassigned raster.
## 4.3 Plot the cropped raster
def pixel_freq_mapbiomas_lulc_raster(geom_name, yr, folder_in_google_drive):
reassigned_raster_path = f'path_to_your_google_drive/{folder_in_google_drive}/cropped_{geom_name}_{yr}.tif'with rasterio.open(reassigned_raster_path) as src:
raster_data = src.read(1)
unique_values = np.unique(raster_data)
total_non_zero = np.sum(raster_data != 255) # Count the entire variety of non-zero pixels
for value in unique_values:
if value != 255: # Exclude no data (= 255)
count = np.sum(raster_data == value) # Count the variety of pixels with the worth
share = count / total_non_zero # Calculate the share of the worth
share = share.round(3)
print(f"Value: {value}, Count: {count}, Share: {share}")
teste3 = pixel_freq_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive)
Lastly, we plot it on a map. You may change the arguments below to regulate traits like colours, labels, legend position, font sizes, etc. Also, there may be an option to decide on the format wherein you wish to save the info (normally PDF or PNG). PDFs are heavier and preserve resolution, while PNGs are lighter but have lower resolution.
## 4.4 Plot the cropped raster
def plot_mapbiomas_lulc_raster(geom_name, yr, folder_in_google_drive,driver):
reassigned_raster_path = f'/Users/vhpf/Library/CloudStorage/GoogleDrive-vh.pires03@gmail.com/My Drive/{folder_in_google_drive}/cropped_{geom_name}_{yr}.tif'
with rasterio.open(reassigned_raster_path) as src:
raster_data = src.read(1)# Define the colours for every class - notice it's essential follow the identical order because the values
values = [3, 10, 15, 18, 22, 26, 255] # Should be the identical of the mapbiomas_categories dictionary
colors_list = ['#6a994e', '#a7c957', '#c32f27', '#dda15e', '#6c757d', '#0077b6','#FFFFFF'] # Set your colours
labels = ['Natural Forest', 'Other Natural Vegetation', 'Pasture', 'Agriculture', 'Others', 'Water', 'No data'] # Set your labels
cmap = colours.ListedColormap(colors_list) # Create a colormap (cmap) with the colours
bounds = values + [256] # Add a worth to the top of the list to incorporate the last color
norm = colours.BoundaryNorm(bounds, cmap.N) # Normalize the colormap to the values
img = plt.imshow(raster_data, interpolation='nearest', cmap=cmap, norm=norm) # Plot the info with the colormap
legend_patches = [mpatches.Patch(color=colors_list[i], label=labels[i]) for i in range(len(values)-1)] # Create the legend patches without the last one (255 = no data)
# Create the legend
plt.legend(handles = legend_patches, # Add the legend patches
bbox_to_anchor = (0.5, -0.02), # Place the legend below the plot
loc = 'upper center', # Place the legend within the upper center
ncol = 3, # Variety of columns
fontsize = 9, # Font size
handlelength=1.5,# Length of the legend handles
frameon=False) # Remove the frame across the legend
plt.axis('off') # Remove the axis
geom_name_title = inflection.titleize(geom_name)
plt.title(f'Land Use in {geom_name_title} ({yr})', fontsize=20) # Add title
saving_path = f'figures/{geom_name}_{yr}.{driver}'
plt.savefig(saving_path, format=driver, dpi=1800) # Reserve it as a .pdf or .png on the figures folder of your project
plt.show()
teste4 = plot_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive, 'png')
Finally, here’s an example of tips on how to use the functions and create a loop to get the LULC evolution for Porto Acre (AC) since 1985. That’s one other city within the AMACRO region with soaring deforestation rates.
## 4.5 Do it in only one function - recall to avoid wasting rasters (4.1) before
def clean_mapbiomas_lulc_raster(geom, geom_name, yr, folder_in_google_drive,driver):
crop_mapbiomas_lulc_raster(geom, geom_name, yr, folder_in_google_drive)
plot_mapbiomas_lulc_raster(geom_name, yr, folder_in_google_drive,driver)
print(f"MapBiomas LULC raster for {geom_name} in {yr} cropped and plotted!")
## 4.6 Run it for multiple geometries for multiple years### 4.6.1 First, save rasters for multiple geometries and years
cities_list = ['Porto Acre'] # Cities to be analyzed - check whether there are two cities in Brazil with the identical name
years = range(1985,2023) # Years to be analyzed (first yr in MapBiomas LULC == 1985)
brazilian_municipalities = gpd.read_file('municipios/file.shp', engine='pyogrio', use_arrow=True) # Read the shapefile - you need to use every other shapefile here
brazilian_municipalities = brazilian_municipalities.clean_names()
brazilian_municipalities.crs = 'EPSG:4326' # Set the CRS to WGS84 (this project default one, change if needed)
selected_cities = brazilian_municipalities.query('nm_mun in @cities_list') # Filter the geometry for the chosen cities
selected_cities = selected_cities.reset_index(drop=True) # Reset the index
cities_ufs = [] # Create list to append the complete names of the cities with their UF (state abbreviation, in portuguese)
nrows = len(selected_cities)
for i in range(nrows):
city = selected_cities.iloc[i]
city_name = city['nm_mun']
city_uf = city['sigla_uf']
cities_ufs.append(f"{city_name} - {city_uf}")
folder_in_google_drive = 'tutorial_mapbiomas_gee' # Folder in Google Drive to avoid wasting the rasters
for city in cities_list:
for yr in years:
city_geom = selected_cities.query(f'nm_mun == "{city}"').geometry.iloc[0] # Get the geometry of the town
geom_name = unidecode.unidecode(city) # Remove latin-1 characters from the town name - GEE doesn`t allow them
get_mapbiomas_lulc_raster(city_geom, geom_name, yr, folder_in_google_drive) # Run the function for every city and yr
### 4.6.2 Second, crop and plot the rasters for multiple geometries and years - Be certain that you may have enough space in your Google Drive and all rasters are there
for city in cities_list:
for yr in years:
city_geom = selected_cities.query(f'nm_mun == "{city}"').geometry.iloc[0]
geom_name = unidecode.unidecode(city)
clean_mapbiomas_lulc_raster(city_geom, geom_name, yr, folder_in_google_drive,'png') # Run the function for every city and yr
gc.collect()
We are going to finish the tutorial by making a short video showing the evolution of deforestation within the municipality during the last 4 many years. Note that you could extend the evaluation to multiple cities and choose specific years for the evaluation. Be happy to customize the algorithm as needed.
## 4.7 Make a clip with LULC evolution
img_folder = 'figures/porto_acre_lulc' # I created a folder to avoid wasting the photographs of the LULC evolution for Porto Acre inside project_path/figures
img_files = sorted([os.path.join(img_folder, f) for f in os.listdir(img_folder) if f.endswith('.png')]) # Gets all the photographs within the folder that end with .png - be certain that you simply have the specified images within the folderclip = ImageSequenceClip(img_files, fps=2) # 2 FPS, 0.5 second between frames
output_file = 'figures/clips/porto_acre_lulc.mp4' # Save clip on the clips folder
clip.write_videofile(output_file, codec='libx264') # It takes some time to create the video (3m30s in my pc)
