Python + Google Earth Engine. The best way to clear MapBiomas LULC rasters for… | by Vinícius Hector | Jul, 2024
# 1. Mission Setup
Initially, we have to load libraries. Be sure that all of them are correctly put in. Additionally, I’m utilizing Python 3.12.3.
## 1.1 Load libraries# If it's worthwhile to set up any library, run under:
# pip set up library1 library2 library3 ...
# Fundamental 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 knowledge 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 Information Manipulation and Visualization
import rasterio # For raster manipulation
from rasterio.masks import masks # For raster knowledge manipulation
from rasterio.plot import present # For raster knowledge visualization
# Plotting and Visualization
import matplotlib.pyplot as plt # For plotting and knowledge 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
# Information Storage and Manipulation
import pyarrow # For environment friendly knowledge storage and manipulation
# Video Making
from moviepy.editor import ImageSequenceClip # For creating movies (part 4.7 solely) - examine this you probably have points: https://github.com/kkroening/ffmpeg-python
Then, be sure you have a folder for this undertaking. All assets and outputs will likely be saved there. This folder will be positioned in your native drive, a cloud-based storage resolution, or in a selected folder on Google Drive the place you’ll save the rasters retrieved utilizing the GEE API.
When operating your code, make certain to vary the handle under to your undertaking path. Home windows customers, at all times 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 should be in
os.chdir(project_path) # All assets on the undertaking 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 perform 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 perform 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 instance to create instinct of the method, we are going to 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 usually 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 could be a sluggish course of, I like to recommend utilizing it provided that it’s worthwhile to take care of just a few or small areas for just a few years. Massive areas might take hours to save lots of on Google Drive, so I like to recommend downloading the heavy LULC information for the entire nation after which cleansing them, as we are going to do in a future put up.
To run the code, first obtain and save IBGE’s¹ Brazilian cities shapefiles (choose Brazil > Municipalities). Bear in mind, you should use any shapefile in Brazil to carry out this algorithm.
## 2.1 Get the geometry of the realm of curiosity (Acrelândia, AC)
brazilian_municipalities = gpd.read_file('municipios/file.shp', engine='pyogrio', use_arrow=True) # Learn the shapefile - you should use every other shapefile right here. Shapefiles should be in your undertaking folder, as set in 1.2
brazilian_municipalities = brazilian_municipalities.clean_names() # Clear the column names (take away particular characters, areas, and many others.)
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 (will be every other 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 now have the shapefile we need to examine correctly saved, we are going to create a bounding field round it to crop the MapBiomas full raster. Then, we are going to put it aside the GEE Python API.
## 2.3 Set the bounding field (bbox) for the realm 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 are going to entry the MapBiomas Google Earth Engine API. First, we have to create a cloud project on GEE utilizing a Google Account. Ensure you have sufficient house in your Google Drive account.
Then, we have to authenticate the GEE Python API (solely as soon as). If you’re a VSCode person, 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 could barely modify this script to work with different property 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 = 'initiatives/mapbiomas-workspace/public/collection8/mapbiomas_collection80_integration_v1'
asset_properties = ee.knowledge.getAsset(mapbiomas_asset) # Examine the asset's properties
asset_properties # See properties
Right here, every band represents the LULC knowledge for a given yr. Guarantee that 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 operating.
I attempted to reserve it domestically, however it appears like it’s worthwhile to save GEE rasters at Google Drive (let me know if you understand how to do it domestically). That is probably the most time-consuming a part of the code. For giant ROIs, this may take hours. Examine your GEE process supervisor to see if the rasters had been correctly loaded to the vacation spot folder.
## 2.6 Export it to your Google Drive (guarantee you may have house there and that it's correctly arrange)# Obs 1: Recall it's worthwhile to authenticate the session, because it was accomplished on 2.4
# Obs 2: Guarantee you may have sufficient house on Google Drive. Free model solely provides 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', # Activity description
folder='tutorial_mapbiomas_gee', # Change this to the folder in your Google Drive the place you need 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 process
export_task.begin()
# 3. Plot the Map
Now we now have a raster with LULC knowledge for a bounding field round Acrelândia in 2022. That is saved on the handle under (at Google Drive). First, let’s see the way it appears.
## 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 knowledge
with rasterio.open(file_path) as src:
knowledge = src.learn(1)
print(src.meta)
print(src.crs)
present(knowledge)
In MapBiomas LULC Assortment 8, every pixel represents a selected land use kind in line with this list. As an example, ‘3’ means ‘Pure Forest’, ‘15’ means ‘Pasture’, and ‘0’ means ‘No knowledge’ (pixels within the raster not throughout the Brazilian borders).
We’ll discover the info we now have earlier than plotting it.
## 3.2 See distinctive values
unique_values = np.distinctive(knowledge)
unique_values # Returns distinctive pixels values within the raster# 0 = no knowledge, elements of the picture outdoors Brazil
## 3.3 See the frequency of every class (besides 0 - no knowledge)
unique_values, counts = np.distinctive(knowledge[data != 0], return_counts=True) # Get the distinctive values and their counts (besides zero)
pixel_counts = pd.DataFrame({'worth': unique_values, 'depend': 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 component represents how every tiny 30m x 30m piece of land is used.
## 3.4 See the precise raster (a matrix through which every component represents a pixel worth - land use code on this case)
knowledge
Now, we have to manage our raster knowledge. As an alternative of categorizing every pixel by actual land use, we are going to categorize them extra broadly. We’ll divide pixels into pure forest, pure non-forest vegetation, water, pasture, agriculture, and different makes use of. Particularly, we’re enthusiastic about monitoring the conversion of pure forest to pasture. To realize this, we are going to reassign pixel values primarily based 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 outdated raster was saved on Google Drive (after being downloaded utilizing GEE’s API), whereas the brand new one will likely be saved within the undertaking folder (in my case, a OneDrive folder on my PC, as set in part 1.2). From right here onwards, we are going to use solely the reassigned raster to plot the info.
That is the primary a part of the script. When you’ve got doubts about what is going on right here (cropping for Acrelândia after which reassigning pixels to broader classes), I like to recommend operating 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 will likely 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 will likely 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 will likely be reassigned to 18 (Agriculture)
# Water ( = 26)
26:26, 33:26, 31:26, # That's, values 26, 33, and 31 will likely 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 will likely be reassigned to 22 (Different)
# No knowledge (= 255)
0:255 # That's, values 0 will likely be reassigned to 255 (No knowledge)
}
## 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 undertaking 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 the town limits (will likely be assigned to no knowledge = 255)
out_image, out_transform = rasterio.masks.masks(src, [city_geom], crop=True)
out_meta.replace({
"top": out_image.form[1],
"width": out_image.form[2],
"rework": 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 stuffed with zeros to be modified
# 3.5.2. Reassign every pixel primarily based on the mapbiomas_categories dictionary
for original_value, new_value in mapbiomas_categories.gadgets():
masks = (raster_array == original_value) # Create a boolean masks for the unique worth (True = Change, False = Do not change)
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, depend=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) # Depend the full variety of non-zero pixelsfor worth in unique_values:
if worth != 255: # Exclude no knowledge (= 255)
depend = np.sum(raster_data == worth) # Depend the variety of pixels with the worth
share = depend / total_non_zero # Calculate the share of the worth
share = share.spherical(3) # Spherical to three decimal locations
print(f"Worth: {worth}, Depend: {depend}, Share: {share}")
Now we plot the info with generic colours. We’ll improve the map later, however that is only a first (or second?) look. Discover that I particularly set 255 (= no knowledge, 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:
knowledge = src.learn(1) # Learn the raster knowledge
unique_values = np.distinctive(knowledge) # Get the distinctive values within the rasterplt.determine(figsize=(10, 8)) # Set the determine dimension
cmap = plt.cm.viridis # Utilizing Viridis colormap
norm = Normalize(vmin=knowledge.min(), vmax=26) # Normalize the info to the vary of the colormap (max = 26, water)
masked_data = np.ma.masked_where(knowledge == 255, knowledge) # Masks no knowledge values (255)
plt.imshow(masked_data, cmap=cmap, norm=norm) # Plot the info with the colormap
plt.colorbar(label='Worth') # Add a colorbar with the values
plt.present()
Now it’s time to create a good looking map. I’ve chosen Matplotlib as a result of I need static maps. In the event you desire interactive choropleths, you should use Plotly.
For additional particulars on choropleths with Matplotlib, examine its documentation, GeoPandas guide, and the nice Yan Holtz’s Python Graph Gallery — the place I get a lot of the inspiration and instruments for DataViz in Python. Additionally, for lovely shade palettes, coolors.co is a wonderful useful resource.
Ensure you have all knowledge visualization libraries correctly loaded to run the code under. I additionally tried to vary the order of patches, however I didn’t know . Let me know when you learn how to do it.
## 3.8 Plot the reassigned raster with customized colours# Outline the colours for every class - discover it's worthwhile to observe the identical order because the values and should be numerically rising or lowering (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 worth to the tip of the record 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 info 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 knowledge)
# 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 dimension
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) # Put it aside as a PDF on the figures folder
plt.present()
4. Standardized Capabilities
Now that we now have constructed instinct on obtain, save, clear, and plot MapBiomas LULC rasters. It’s time to generalize the method.
On this part, we are going to outline capabilities to automate these steps for any form and any yr. Then, we are going to execute these capabilities in a loop to research a selected metropolis — Porto Acre, AC — from 1985 to 2022. Lastly, we are going to make a video illustrating the LULC evolution in that space over the required 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 perform 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 = 'initiatives/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}", # Activity description
folder=folder_in_google_drive, # Change this to the folder in your Google Drive the place you need 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 might 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 undertaking 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’ll be able to change `reassigned_raster_path` to reserve it anyplace else. In the event you change it, make certain to replace the remainder of the code accordingly.
Additionally, you’ll be able to reassign pixels as wanted by modifying the mapbiomas_categories dict. The left digit represents the unique pixel values, and the proper 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 knowledge (= 255)
0:255
} # You possibly can 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 the town limits (will likely be assigned to no knowledge = 0)
out_image, out_transform = rasterio.masks.masks(src, [my_geom], crop=True)
out_meta.replace({
"top": out_image.form[1],
"width": out_image.form[2],
"rework": 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 stuffed with zeros to be modified
# Reassign every pixel primarily based on the mapbiomas_categories dictionary
for original_value, new_value in mapbiomas_categories.gadgets():
masks = (raster_array == original_value) # Create a boolean masks for the unique worth (True = Change, False = Do not change)
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, depend=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) # Depend the full variety of non-zero pixels
for worth in unique_values:
if worth != 255: # Exclude no knowledge (= 255)
depend = np.sum(raster_data == worth) # Depend the variety of pixels with the worth
share = depend / total_non_zero # Calculate the share of the worth
share = share.spherical(3)
print(f"Worth: {worth}, Depend: {depend}, Share: {share}")
teste3 = pixel_freq_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive)
Lastly, we plot it on a map. You possibly can change the arguments under to regulate traits like colours, labels, legend place, font sizes, and many others. Additionally, there’s an possibility to decide on the format through which you need to save the info (normally 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 it's worthwhile to 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 worth to the tip of the record 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 info 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 knowledge)
# 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 dimension
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) # Put it aside as a .pdf or .png on the figures folder of your undertaking
plt.present()
teste4 = plot_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive, 'png')
Lastly, right here’s an instance of 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 perform - recall to save lots of 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 - examine 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 should use every other shapefile right here
brazilian_municipalities = brazilian_municipalities.clean_names()
brazilian_municipalities.crs = 'EPSG:4326' # Set the CRS to WGS84 (this undertaking 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 record 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 save lots of 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 the town
geom_name = unidecode.unidecode(metropolis) # Take away latin-1 characters from the town title - GEE doesn`t permit them
get_mapbiomas_lulc_raster(city_geom, geom_name, yr, folder_in_google_drive) # Run the perform 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 house 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 perform for every metropolis and yr
gc.gather()
We’ll end the tutorial by creating a brief video displaying the evolution of deforestation within the municipality over the past 4 a long time. Notice that you could prolong 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 save lots of 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 you 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 laptop)
