# 1. Venture Setup
To begin with, we have to load libraries. Be certain all of them are correctly put in. Additionally, I’m utilizing Python 3.12.3.
## 1.1 Load libraries# If you could set up any library, run under:
# pip set up library1 library2 library3 ...
# Primary Libraries
import os # For file operations
import gc # For rubbish assortment, it avoids RAM reminiscence points
import numpy as np # For coping with matrices
import pandas as pd # For coping with dataframes
import janitor # For information cleansing (primarily column names)
import numexpr # For quick pd.question() manipulation
import inflection # For string manipulation
import unidecode # For string manipulation
# Geospatial Libraries
import geopandas as gpd # For coping with shapefiles
import pyogrio # For quick .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, Level, MultiPolygon, LineString # For geometry manipulation
# Raster Knowledge Manipulation and Visualization
import rasterio # For raster manipulation
from rasterio.masks import masks # For raster information manipulation
from rasterio.plot import present # For raster information visualization
# Plotting and Visualization
import matplotlib.pyplot as plt # For plotting and information visualization
from matplotlib.colours import ListedColormap, Normalize # For shade manipulation
import matplotlib.colours as colours # For shade manipulation
import matplotlib.patches as mpatches # For creating patch objects
import matplotlib.cm as cm # For colormaps
# Knowledge Storage and Manipulation
import pyarrow # For environment friendly information storage and manipulation
# Video Making
from moviepy.editor import ImageSequenceClip # For creating movies (part 4.7 solely) - verify this when you have points: https://github.com/kkroening/ffmpeg-python
Then, be sure to have a folder for this mission. All assets and outputs can be saved there. This folder could be positioned in your native drive, a cloud-based storage answer, or in a selected folder on Google Drive the place you’ll save the rasters retrieved utilizing the GEE API.
When working your code, be certain that to alter the tackle under to your mission path. Home windows customers, all the time keep in mind to make use of as an alternative of /.
# 1.2 Set working listing
project_path = 'path_to_your_project_folder' # The place you'll save all outcomes and assets have to be in
os.chdir(project_path) # All assets on the mission are relative to this path# 1.3 Additional settings
pd.set_option('compute.use_numexpr', True) # Use numexpr for quick pd.question() manipulation
Lastly, this operate is beneficial for plotting geometries over OpenStreetMap (OSM). It’s notably useful when working with unknown shapefiles to make sure accuracy and keep away from errors.
## 1.4 Set operate to plot geometries over an OSM
def plot_geometries_on_osm(geometries, zoom_start=10):# Calculate the centroid of all geometries to heart 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 across the common centroid
map = folium.Map(location=[avg_y, avg_x], zoom_start=zoom_start)
# Add every 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 Instance: Acrelândia (AC) in 2022
For example to create instinct of the method, we’ll save, clear, and plot land use in Acrelândia (AC) in 2022. It’s a metropolis in the midst of the AMACRO area (the three-state border of Amazonas, Acre, and Rondônia), the place the customarily untouched forest is being quickly destroyed.
On this part, I’ll clarify step-by-step of the script, after which standardize the method to run it for a number of locations over a number of years. Since saving giant rasters utilizing the API is usually a sluggish course of, I like to recommend utilizing it provided that you could cope with a number of or small areas for a number of years. Massive areas could take hours to avoid wasting on Google Drive, so I like to recommend downloading the heavy LULC recordsdata for the entire nation after which cleansing them, as we’ll do in a future publish.
To run the code, first obtain and save IBGE’s¹ Brazilian cities shapefiles (choose Brazil > Municipalities). Bear in mind, you need to use any shapefile in Brazil to carry out this algorithm.
## 2.1 Get the geometry of the world of curiosity (Acrelândia, AC)
brazilian_municipalities = gpd.read_file('municipios/file.shp', engine='pyogrio', use_arrow=True) # Learn the shapefile - you need to use another shapefile right here. Shapefiles have to be in your mission folder, as set in 1.2
brazilian_municipalities = brazilian_municipalities.clean_names() # Clear the column names (take away particular characters, areas, and so on.)
brazilian_municipalities.crs = 'EPSG:4326' # Set the CRS to WGS84 (MapBiomas makes use of this CRS)
brazilian_municipalities
## 2.2 Get geometry for Acrelândia, AC
metropolis = brazilian_municipalities.question('nm_mun == "Acrelândia"') # Filter the geometry for Acrelândia, AC (could be another metropolis or set of cities)
city_geom = metropolis.geometry.iloc[0] # Get the geometry of Acrelândia, AC
city_geom # See the geometry form
As soon as we’ve the shapefile we wish to research correctly saved, we’ll create a bounding field round it to crop the MapBiomas full raster. Then, we’ll put it aside the GEE Python API.
## 2.3 Set the bounding field (bbox) for the world of curiosity
bbox = city_geom.bounds # Get the bounding field of the geometry
bbox = Polygon([(bbox[0], bbox[1]), (bbox[0], bbox[3]), (bbox[2], bbox[3]), (bbox[2], bbox[1])]) # Convert the bounding field to a Polygonbbox_xmin = bbox.bounds[0] # Get the minimal x coordinate of the bounding field
bbox_ymin = bbox.bounds[1] # Get the minimal y coordinate of the bounding field
bbox_xmax = bbox.bounds[2] # Get the utmost x coordinate of the bounding field
bbox_ymax = bbox.bounds[3] # Get the utmost y coordinate of the bounding field
bbox # See bbox round Acrelândia form
# Plot the bounding field and the geometry of Acrelandia over an OSM map
plot_geometries_on_osm([bbox, city_geom], zoom_start=10)
Now, we’ll entry the MapBiomas Google Earth Engine API. First, we have to create a cloud mission on GEE utilizing a Google Account. Ensure you have sufficient area in your Google Drive account.
Then, we have to authenticate the GEE Python API (solely as soon as). If you’re a VSCode consumer, discover that the token insertion field seems within the high proper nook of the IDE.
All photos from the MapBiomas LULC Assortment can be found in the identical asset. Discover that you would be able to barely modify this script to work with different belongings within the GEE catalog and different MapBiomas collections.
## 2.4 Acess MapBiomas Assortment 8.0 utilizing 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
# Outline the MapBiomas Assortment 8.0 asset ID - retrieved from https://brasil.mapbiomas.org/en/colecoes-mapbiomas/
mapbiomas_asset = 'tasks/mapbiomas-workspace/public/collection8/mapbiomas_collection80_integration_v1'
asset_properties = ee.information.getAsset(mapbiomas_asset) # Verify the asset's properties
asset_properties # See properties
Right here, every band represents the LULC information for a given yr. Make it possible for the code under is correctly written. This selects the picture for the specified yr after which crops the uncooked raster for a bounding field across the area of curiosity (ROI) — Acrelândia, AC.
## 2.5 Filter the gathering for 2022 and crop the gathering to a bbox round 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.Picture(mapbiomas_asset) # Get the pictures of MapBiomas Assortment 8.0
mapbiomas2022 = mapbiomas_image.choose(band_id) # Choose the picture for 2022
roi = ee.Geometry.Rectangle([bbox_xmin, bbox_ymin, bbox_xmax, bbox_ymax]) # Set the Area of Curiosity (ROI) to the bbox round Acrelândia, AC - set in 2.3
image_roi = mapbiomas2022.clip(roi) # Crop the picture to the ROI
Now, we save the cropped raster on Google Drive (in my case, into the ‘tutorial_mapbiomas_gee’ folder). Ensure you have created the vacation spot folder in your drive earlier than working.
I attempted to reserve it regionally, however it appears to be like like you could save GEE rasters at Google Drive (let me know if you know the way to do it regionally). That is essentially the most time-consuming a part of the code. For big ROIs, this would possibly take hours. Verify your GEE activity supervisor to see if the rasters have been correctly loaded to the vacation spot folder.
## 2.6 Export it to your Google Drive (guarantee you've area there and that it's correctly arrange)# Obs 1: Recall you could authenticate the session, because it was finished on 2.4
# Obs 2: Guarantee you've sufficient area on Google Drive. Free model solely offers 15 Gb of storage.
export_task = ee.batch.Export.picture.toDrive(
picture=image_roi, # Picture to export to Google Drive as a GeoTIFF
description='clipped_mapbiomas_collection8_acrelandia_ac_2022', # Job description
folder='tutorial_mapbiomas_gee', # Change this to the folder in your Google Drive the place you wish to save the file
fileNamePrefix='acrelandia_ac_2022', # File title (change it if you wish to)
area=roi.getInfo()['coordinates'], # Area to export the picture
scale=30,
fileFormat='GeoTIFF'
)
# Begin the export activity
export_task.begin()
# 3. Plot the Map
Now we’ve a raster with LULC information for a bounding field round Acrelândia in 2022. That is saved on the tackle under (at Google Drive). First, let’s see the way it appears to be like.
## 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 information
with rasterio.open(file_path) as src:
information = src.learn(1)
print(src.meta)
print(src.crs)
present(information)
In MapBiomas LULC Assortment 8, every pixel represents a selected land use sort in line with this checklist. As an example, ‘3’ means ‘Pure Forest’, ‘15’ means ‘Pasture’, and ‘0’ means ‘No information’ (pixels within the raster not throughout the Brazilian borders).
We are going to discover the information we’ve earlier than plotting it.
## 3.2 See distinctive values
unique_values = np.distinctive(information)
unique_values # Returns distinctive pixels values within the raster# 0 = no information, elements of the picture outdoors Brazil
## 3.3 See the frequency of every class (besides 0 - no information)
unique_values, counts = np.distinctive(information[data != 0], return_counts=True) # Get the distinctive values and their counts (besides zero)
pixel_counts = pd.DataFrame({'worth': unique_values, 'rely': counts})
pixel_counts['share'] = (pixel_counts['count'] / pixel_counts['count'].sum())*100
pixel_counts
On the finish of the day, we’re working with a big matrix through which every ingredient represents how every tiny 30m x 30m piece of land is used.
## 3.4 See the precise raster (a matrix through which every ingredient represents a pixel worth - land use code on this case)
information
Now, we have to set up our raster information. As an alternative of categorizing every pixel by actual land use, we’ll categorize them extra broadly. We are going to divide pixels into pure forest, pure non-forest vegetation, water, pasture, agriculture, and different makes use of. Particularly, we’re excited by monitoring the conversion of pure forest to pasture. To attain this, we’ll reassign pixel values based mostly on the mapbiomas_categories dict under, which follows with MapBiomas’ land use and land cowl (LULC) categorization.
The code under crops the raster to Acrelândia’s limits and reassigns pixels in line with the mapbiomas_categories dict. Then, it saves it as a brand new raster at ‘reassigned_raster_path’. Notice that the previous raster was saved on Google Drive (after being downloaded utilizing GEE’s API), whereas the brand new one can be saved within the mission folder (in my case, a OneDrive folder on my PC, as set in part 1.2). From right here onwards, we’ll use solely the reassigned raster to plot the information.
That is the primary a part of the script. In case you have doubts about what is going on right here (cropping for Acrelândia after which reassigning pixels to broader classes), I like to recommend working it and printing outcomes 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 can be reassigned to three (Forest)
# Different Non-Forest Pure 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 can be reassigned to 10 (Different Non-Forest Pure 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 can be reassigned to 18 (Agriculture)
# Water ( = 26)
26:26, 33:26, 31:26, # That's, values 26, 33, and 31 can be reassigned to 26 (Water)
# Different (= 22)
22:22, 23:22, 24:22, 30:22, 25:22, 27:22, # That's, values 22, 23, 24, 30, 25, and 27 can be reassigned to 22 (Different)
# No information (= 255)
0:255 # That's, values 0 can be reassigned to 255 (No information)
}
## 3.5 Reassing pixels values to the MapBiomas customized basic classes 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' # Someplace within the mission folder set at 1.2with rasterio.open(original_raster_path) as src:
raster_array = src.learn(1)
out_meta = src.meta.copy() # Get metadata from the unique raster
# 3.5.1. Crop (or masks) the raster to the geometry of city_geom (on this case, Acrelandia, AC) and thus take away pixels outdoors town limits (can be assigned to no information = 255)
out_image, out_transform = rasterio.masks.masks(src, [city_geom], crop=True)
out_meta.replace({
"peak": out_image.form[1],
"width": out_image.form[2],
"remodel": out_transform
}) # Replace 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 every pixel based mostly on the mapbiomas_categories dictionary
for original_value, new_value in mapbiomas_categories.objects():
masks = (raster_array == original_value) # Create a boolean masks for the unique worth (True = Change, False = Do not exchange)
modified_raster[mask] = new_value # Change the unique values with the brand new values, when wanted (that's, when the masks is True)
out_meta = src.meta.copy() # Get metadata from the unique raster
out_meta.replace(dtype=rasterio.uint8, rely=1) # Replace 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.learn(1)
unique_values = np.distinctive(raster_data)
total_non_zero = np.sum(raster_data != 255) # Rely the entire variety of non-zero pixelsfor worth in unique_values:
if worth != 255: # Exclude no information (= 255)
rely = np.sum(raster_data == worth) # Rely the variety of pixels with the worth
share = rely / total_non_zero # Calculate the share of the worth
share = share.spherical(3) # Spherical to three decimal locations
print(f"Worth: {worth}, Rely: {rely}, Share: {share}")
Now we plot the information with generic colours. We are going to improve the map later, however that is only a first (or second?) look. Discover that I particularly set 255 (= no information, pixels outdoors 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:
information = src.learn(1) # Learn the raster information
unique_values = np.distinctive(information) # Get the distinctive values within the rasterplt.determine(figsize=(10, 8)) # Set the determine measurement
cmap = plt.cm.viridis # Utilizing Viridis colormap
norm = Normalize(vmin=information.min(), vmax=26) # Normalize the information to the vary of the colormap (max = 26, water)
masked_data = np.ma.masked_where(information == 255, information) # Masks no information values (255)
plt.imshow(masked_data, cmap=cmap, norm=norm) # Plot the information with the colormap
plt.colorbar(label='Worth') # Add a colorbar with the values
plt.present()
Now it’s time to create a wonderful map. I’ve chosen Matplotlib as a result of I would like static maps. Should you choose interactive choropleths, you need to use Plotly.
For additional particulars on choropleths with Matplotlib, verify its documentation, GeoPandas information, and the good Yan Holtz’s Python Graph Gallery — the place I get a lot of the inspiration and instruments for DataViz in Python. Additionally, for stunning shade palettes, coolors.co is a wonderful useful resource.
Ensure you have all information visualization libraries correctly loaded to run the code under. I additionally tried to alter the order of patches, however I didn’t know methods to. Let me know in case you learn how to do it.
## 3.8 Plot the reassigned raster with customized colours# Outline the colours for every class - discover you could observe the identical order because the values and have to be numerically growing or reducing (nonetheless must learn how to unravel 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 price to the tip of the checklist to incorporate the final shade
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 information with the colormap
legend_patches = [mpatches.Patch(color=colors_list[i], label=labels[i]) for i in vary(len(values)-1)] # Create the legend patches withou the final one (255 = no information)
# Create the legend
plt.legend(handles = legend_patches, # Add the legend patches
bbox_to_anchor = (0.5, -0.02), # Place the legend under the plot
loc = 'higher heart', # Place the legend within the higher heart
ncol = 3, # Variety of columns
fontsize = 9, # Font measurement
handlelength=1,# Size of the legend handles
frameon=False) # Take away the body across the legend
plt.axis('off') # Take away 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.present()
4. Standardized Features
Now that we’ve constructed instinct on methods to obtain, save, clear, and plot MapBiomas LULC rasters. It’s time to generalize the method.
On this part, we’ll outline capabilities to automate these steps for any form and any yr. Then, we’ll execute these capabilities in a loop to investigate a selected metropolis — Porto Acre, AC — from 1985 to 2022. Lastly, we’ll make a video illustrating the LULC evolution in that space over the desired interval.
First, save a bounding field (bbox) across the Area of Curiosity (ROI). You solely must enter the specified geometry and specify the yr. This operate 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 field of the geometry
bbox = Polygon([(bbox[0], bbox[1]), (bbox[0], bbox[3]), (bbox[2], bbox[3]), (bbox[2], bbox[1])]) # Convert the bounding field to a Polygon
bbox_xmin = bbox.bounds[0] # Get the minimal x coordinate of the bounding field
bbox_ymin = bbox.bounds[1] # Get the minimal y coordinate of the bounding field
bbox_xmax = bbox.bounds[2] # Get the utmost x coordinate of the bounding field
bbox_ymax = bbox.bounds[3] # Get the utmost y coordinate of the bounding field
mapbiomas_asset = 'tasks/mapbiomas-workspace/public/collection8/mapbiomas_collection80_integration_v1'
band_id = f'classification_{yr}'
mapbiomas_image = ee.Picture(mapbiomas_asset) # Get the pictures of MapBiomas Assortment 8.0
mapbiomas_data = mapbiomas_image.choose(band_id) # Choose the picture for 2022
roi = ee.Geometry.Rectangle([bbox_xmin, bbox_ymin, bbox_xmax, bbox_ymax]) # Set the Area of Curiosity (ROI) to the bbox across the desired geometry
image_roi = mapbiomas_data.clip(roi) # Crop the picture to the ROI
export_task = ee.batch.Export.picture.toDrive(
picture=image_roi, # Picture to export to Google Drive as a GeoTIFF
description=f"save_bbox_around_{geom_name}_in_{yr}", # Job description
folder=folder_in_google_drive, # Change this to the folder in your Google Drive the place you wish to save the file
fileNamePrefix=f"{geom_name}_{yr}", # File title
area=roi.getInfo()['coordinates'], # Area to export the picture
scale=30,
fileFormat='GeoTIFF'
)
export_task.begin() # Discover that importing these rasters to Google Drive could take some time, specifically for big areas
# Check it utilizing Rio de Janeiro in 2022
folder_in_google_drive = 'tutorial_mapbiomas_gee'
rio_de_janeiro = brazilian_municipalities.question('nm_mun == "Rio de Janeiro"')
rio_de_janeiro.crs = 'EPSG:4326' # Set the CRS to WGS84 (this mission default one, change if wanted)
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 solely the pixels throughout the geometry and put it aside as a brand new raster.
I selected to reserve it on Google Drive, however you possibly can change `reassigned_raster_path` to reserve it wherever else. Should you change it, be certain that to replace the remainder of the code accordingly.
Additionally, you possibly can reassign pixels as wanted by modifying the mapbiomas_categories dict. The left digit represents the unique pixel values, and the suitable one represents the reassigned (new) 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,
# Different Non-Forest Pure 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,
# Different (= 22)
22:22, 23:22, 24:22, 30:22, 25:22, 27:22,
# No information (= 255)
0:255
} # You'll be able to change this to no matter categorization you need, however simply keep in mind to adapt the colours and labels within the plot
with rasterio.open(original_raster_path) as src:
raster_array = src.learn(1)
out_meta = src.meta.copy() # Get metadata from the unique raster
# Crop the raster to the geometry of my_geom and thus take away pixels outdoors town limits (can be assigned to no information = 0)
out_image, out_transform = rasterio.masks.masks(src, [my_geom], crop=True)
out_meta.replace({
"peak": out_image.form[1],
"width": out_image.form[2],
"remodel": out_transform
}) # Replace 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 every pixel based mostly on the mapbiomas_categories dictionary
for original_value, new_value in mapbiomas_categories.objects():
masks = (raster_array == original_value) # Create a boolean masks for the unique worth (True = Change, False = Do not exchange)
modified_raster[mask] = new_value # Change the unique values with the brand new values, when wanted (that's, when the masks is True)
out_meta = src.meta.copy() # Get metadata from the unique raster
out_meta.replace(dtype=rasterio.uint8, rely=1) # Replace 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.learn(1)
unique_values = np.distinctive(raster_data)
total_non_zero = np.sum(raster_data != 255) # Rely the entire variety of non-zero pixels
for worth in unique_values:
if worth != 255: # Exclude no information (= 255)
rely = np.sum(raster_data == worth) # Rely the variety of pixels with the worth
share = rely / total_non_zero # Calculate the share of the worth
share = share.spherical(3)
print(f"Worth: {worth}, Rely: {rely}, Share: {share}")
teste3 = pixel_freq_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive)
Lastly, we plot it on a map. You’ll be able to change the arguments under to regulate traits like colours, labels, legend place, font sizes, and so on. Additionally, there’s an choice to decide on the format through which you wish to save the information (often PDF or PNG). PDFs are heavier and protect decision, whereas PNGs are lighter however have decrease decision.
## 4.4 Plot the cropped raster
def plot_mapbiomas_lulc_raster(geom_name, yr, folder_in_google_drive,driver):
reassigned_raster_path = f'/Customers/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.learn(1)# Outline the colours for every class - discover you could observe the identical order because the values
values = [3, 10, 15, 18, 22, 26, 255] # Have to 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 price to the tip of the checklist to incorporate the final shade
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 information with the colormap
legend_patches = [mpatches.Patch(color=colors_list[i], label=labels[i]) for i in vary(len(values)-1)] # Create the legend patches with out the final one (255 = no information)
# Create the legend
plt.legend(handles = legend_patches, # Add the legend patches
bbox_to_anchor = (0.5, -0.02), # Place the legend under the plot
loc = 'higher heart', # Place the legend within the higher heart
ncol = 3, # Variety of columns
fontsize = 9, # Font measurement
handlelength=1.5,# Size of the legend handles
frameon=False) # Take away the body across the legend
plt.axis('off') # Take away 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 mission
plt.present()
teste4 = plot_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive, 'png')
Lastly, right here’s an instance of methods to use the capabilities and create a loop to get the LULC evolution for Porto Acre (AC) since 1985. That’s one other metropolis within the AMACRO area with hovering deforestation charges.
## 4.5 Do it in only one operate - recall to avoid wasting rasters (4.1) earlier than
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 a number of geometries for a number of years### 4.6.1 First, save rasters for a number of geometries and years
cities_list = ['Porto Acre'] # Cities to be analyzed - verify whether or not there are two cities in Brazil with the identical title
years = vary(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) # Learn the shapefile - you need to use another shapefile right here
brazilian_municipalities = brazilian_municipalities.clean_names()
brazilian_municipalities.crs = 'EPSG:4326' # Set the CRS to WGS84 (this mission default one, change if wanted)
selected_cities = brazilian_municipalities.question('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 checklist to append the complete names of the cities with their UF (state abbreviation, in portuguese)
nrows = len(selected_cities)
for i in vary(nrows):
metropolis = selected_cities.iloc[i]
city_name = metropolis['nm_mun']
city_uf = metropolis['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 metropolis in cities_list:
for yr in years:
city_geom = selected_cities.question(f'nm_mun == "{metropolis}"').geometry.iloc[0] # Get the geometry of town
geom_name = unidecode.unidecode(metropolis) # Take away latin-1 characters from town title - GEE doesn`t permit them
get_mapbiomas_lulc_raster(city_geom, geom_name, yr, folder_in_google_drive) # Run the operate for every metropolis and yr
### 4.6.2 Second, crop and plot the rasters for a number of geometries and years - Ensure you have sufficient area in your Google Drive and all rasters are there
for metropolis in cities_list:
for yr in years:
city_geom = selected_cities.question(f'nm_mun == "{metropolis}"').geometry.iloc[0]
geom_name = unidecode.unidecode(metropolis)
clean_mapbiomas_lulc_raster(city_geom, geom_name, yr, folder_in_google_drive,'png') # Run the operate for every metropolis and yr
gc.gather()
We are going to end the tutorial by creating a brief video exhibiting the evolution of deforestation within the municipality over the past 4 a long time. Notice that you would be able to lengthen the evaluation to a number of cities and choose particular years for the evaluation. Be happy to customise the algorithm as wanted.
## 4.7 Make a clip with LULC evolution
img_folder = 'figures/porto_acre_lulc' # I created a folder to avoid wasting the pictures 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')]) # Will get all the pictures within the folder that finish with .png - be sure to solely have the specified photos 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 computer)
