[ad_1]
# 1. Thiết lập dự án
Trước hết, chúng ta cần tải thư viện. Đảm bảo tất cả chúng được cài đặt đúng cách. Ngoài ra, tôi đang sử dụng Python 3.12.3.
## 1.1 Load libraries# If you must set up any library, run beneath:
# 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) - verify this when you have points: https://github.com/kkroening/ffmpeg-python
Sau đó, hãy đảm bảo bạn có một thư mục cho dự án này. Tất cả các tài nguyên và đầu ra sẽ được lưu ở đó. Thư mục này có thể nằm trên ổ đĩa cục bộ của bạn, giải pháp lưu trữ dựa trên đám mây hoặc trong một thư mục cụ thể trên Google Drive, nơi bạn sẽ lưu các raster được lấy bằng API GEE.
Khi chạy mã của bạn, hãy đảm bảo thay đổi địa chỉ bên dưới thành đường dẫn dự án của bạn. Người dùng Home windows, hãy luôn nhớ sử dụng thay vì /.
# 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 venture are relative to this path# 1.3 Additional settings
pd.set_option('compute.use_numexpr', True) # Use numexpr for quick pd.question() manipulation
Cuối cùng, chức năng này hữu ích khi vẽ hình học trên OpenStreetMap (OSM). Nó đặc biệt hữu ích khi làm việc với các shapefile không xác định để đảm bảo độ chính xác và tránh sai sót.
## 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 middle 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. Ví dụ đơn lẻ: Acrelândia (AC) năm 2022
Để tạo ra trực giác về quy trình, chúng ta sẽ tiết kiệm, dọn dẹp và lập biểu đồ sử dụng đất trong Acrelândia (AC) vào năm 2022. Đó là một thành phố ở giữa AMACRO khu vực (biên giới ba tiểu bang của Làazonas, AClại, và Rondônia), nơi khu rừng thường không bị tác động đang bị phá hủy nhanh chóng.
Trong phần này, tôi sẽ giải thích từng bước của tập lệnh, sau đó chuẩn hóa quy trình để chạy nó cho nhiều nơi trong nhiều năm. Vì việc lưu các raster lớn bằng API có thể là một quá trình chậm, tôi khuyên bạn chỉ nên sử dụng nó nếu bạn cần xử lý một vài hoặc một số khu vực nhỏ trong vài năm. Các khu vực lớn có thể mất nhiều giờ để lưu trên Google Drive, vì vậy tôi khuyên bạn nên tải xuống các tệp LULC nặng cho toàn bộ quốc gia rồi dọn dẹp chúng, như chúng ta sẽ làm trong một bài đăng trong tương lai.
Để chạy mã, trước tiên hãy tải xuống và lưu IBGE¹ Các thành phố của Brazil shapefile (chọn Brazil > Thành phố). Hãy nhớ rằng bạn có thể sử dụng bất kỳ shapefile nào ở Brazil để thực hiện thuật toán này.
## 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 should use another shapefile right here. Shapefiles should be in your venture folder, as set in 1.2
brazilian_municipalities = brazilian_municipalities.clean_names() # Clear the column names (take away particular characters, areas, and so forth.)
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 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
Sau khi lưu shapefile mà chúng ta muốn nghiên cứu đúng cách, chúng ta sẽ tạo một hộp giới hạn xung quanh nó để cắt toàn bộ raster MapBiomas. Sau đó, chúng ta sẽ lưu nó trong 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)
Bây giờ, chúng ta sẽ truy cập API Google Earth Engine của MapBiomas. Đầu tiên, chúng ta cần tạo một dự án đám mây trên GEE sử dụng Tài khoản Google. Đảm bảo bạn có đủ dung lượng trên tài khoản Google Drive của mình.
Sau đó, chúng ta cần xác thực API Python của GEE (chỉ một lần). Nếu bạn là người dùng VSCode, hãy lưu ý rằng hộp chèn mã thông báo xuất hiện ở góc trên bên phải của IDE.
Tất cả hình ảnh từ Bộ sưu tập LULC của MapBiomas đều có sẵn trong cùng một tài sản. Lưu ý rằng bạn có thể sửa đổi một chút tập lệnh này để làm việc với các tài sản khác trong Danh mục GEE và các bộ sưu tập MapBiomas khác.
## 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
Ở đây, mỗi băng tần biểu diễn dữ liệu LULC cho một năm nhất định. Đảm bảo rằng mã bên dưới được viết đúng. Mã này sẽ chọn hình ảnh cho năm mong muốn và sau đó cắt raster thô để tạo hộp giới hạn xung quanh vùng quan tâm (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
Bây giờ, chúng ta lưu raster đã cắt trên Google Drive (trong trường hợp của tôi, vào thư mục ‘tutorial_mapbiomas_gee’). Đảm bảo bạn đã tạo thư mục đích trong ổ đĩa của mình trước khi chạy.
Tôi đã thử lưu cục bộ, nhưng có vẻ như bạn cần lưu các raster GEE tại Google Drive (hãy cho tôi biết nếu bạn biết cách thực hiện tại địa phương). Đây là phần tốn nhiều thời gian nhất của mã. Đối với các ROI lớn, việc này có thể mất nhiều giờ. Kiểm tra trình quản lý tác vụ GEE của bạn để xem các raster đã được tải đúng cách vào thư mục đích chưa.
## 2.6 Export it to your Google Drive (guarantee you could have house there and that it's correctly arrange)# Obs 1: Recall you must authenticate the session, because it was completed on 2.4
# Obs 2: Guarantee you could 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 wish to save the file
fileNamePrefix='acrelandia_ac_2022', # File identify (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. Vẽ bản đồ
Bây giờ chúng ta có một raster với dữ liệu LULC cho một hộp giới hạn xung quanh Acrelândia vào năm 2022. Dữ liệu này được lưu tại địa chỉ bên dưới (tại Google Drive). Trước tiên, hãy xem nó trông như thế nào.
## 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)
Trong Bộ sưu tập LULC 8 của MapBiomas, mỗi pixel biểu thị một loại sử dụng đất cụ thể theo danh sách này. Ví dụ, ‘3’ có nghĩa là ‘Rừng tự nhiên’, ’15’ có nghĩa là ‘Đồng cỏ’ và ‘0’ có nghĩa là ‘Không có dữ liệu’ (điểm ảnh trong raster không nằm trong biên giới Brazil).
Chúng ta sẽ khám phá dữ liệu mình có trước khi vẽ biểu đồ.
## 3.2 See distinctive values
unique_values = np.distinctive(knowledge)
unique_values # Returns distinctive pixels values within the raster# 0 = no knowledge, components of the picture exterior Brazil
## 3.3 See the frequency of every class (besides 0 - no knowledge)
unique_values, counts = np.distinctive(knowledge(knowledge != 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('rely') / pixel_counts('rely').sum())*100
pixel_counts
Vào cuối ngày, chúng ta sẽ làm việc với một ma trận lớn trong đó mỗi phần tử biểu thị cách sử dụng từng mảnh đất nhỏ có kích thước 30m x 30m.
## 3.4 See the precise raster (a matrix during which every aspect represents a pixel worth - land use code on this case)
knowledge
Bây giờ, chúng ta cần sắp xếp dữ liệu raster của mình. Thay vì phân loại từng pixel theo mục đích sử dụng đất chính xác, chúng ta sẽ phân loại chúng rộng hơn. Chúng ta sẽ chia pixel thành rừng tự nhiên, thảm thực vật tự nhiên không phải rừng, nước, đồng cỏ, nông nghiệpVà khác sử dụng. Cụ thể, chúng tôi quan tâm đến việc theo dõi quá trình chuyển đổi rừng tự nhiên thành đồng cỏ. Để đạt được điều này, chúng tôi sẽ chỉ định lại các giá trị pixel dựa trên mapbiomas_categories dict bên dưới, theo sau là phân loại sử dụng đất và lớp phủ đất (LULC) của MapBiomas.
Mã bên dưới cắt raster theo giới hạn của Acrelândia và chỉ định lại các pixel theo mapbiomas_categories dict. Sau đó, nó lưu nó dưới dạng raster mới tại ‘reassigned_raster_path’. Lưu ý rằng raster cũ đã được lưu trên Google Drive (sau khi được tải xuống bằng API của GEE), trong khi raster mới sẽ được lưu trong thư mục dự án (trong trường hợp của tôi, là thư mục OneDrive trên PC của tôi, như đã đặt trong phần 1.2). Từ đây trở đi, chúng ta sẽ chỉ sử dụng raster được gán lại để vẽ dữ liệu.
Đây là phần chính của kịch bản. Nếu bạn có thắc mắc về những gì đang diễn ra ở đây (cắt xén cho Acrelândia rồi gán lại pixel cho các danh mục rộng hơn), tôi khuyên bạn nên chạy chương trình và in kết quả cho từng bước.
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 venture 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 exterior 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({
"peak": 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 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 = Exchange, False = Do not substitute)
modified_raster(masks) = new_value # Exchange 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) # Depend the full variety of non-zero pixelsfor worth in unique_values:
if worth != 255: # Exclude no knowledge (= 255)
rely = np.sum(raster_data == worth) # Depend 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}, Depend: {rely}, Share: {share}")
Bây giờ chúng ta vẽ dữ liệu bằng màu chung. Chúng ta sẽ cải thiện bản đồ sau, nhưng đây chỉ là cái nhìn đầu tiên (hoặc thứ hai?). Lưu ý rằng tôi đặc biệt đặt 255 (= không có dữ liệu, pixel bên ngoài Acrelândia) thành màu trắng để trực quan hóa tốt hơn.
## 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 information 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 information with the colormap
plt.colorbar(label='Worth') # Add a colorbar with the values
plt.present()
Bây giờ là lúc tạo một bản đồ đẹp. Tôi đã chọn Matplotlib vì tôi muốn có bản đồ tĩnh. Nếu bạn thích bản đồ choropleth tương tác, bạn có thể sử dụng Cốt truyện.
Để biết thêm chi tiết về choropleths với Matplotlib, hãy kiểm tra tài liệu của nó, Hướng dẫn GeoPandasvà Yan Holtz vĩ đại Thư viện đồ thị Python — nơi tôi có được nhiều cảm hứng và công cụ cho DataViz trong Python. Ngoài ra, đối với bảng màu đẹp, coolors.co là một nguồn tài nguyên tuyệt vời.
Hãy đảm bảo bạn đã tải đúng tất cả các thư viện trực quan hóa dữ liệu để chạy mã bên dưới. Tôi cũng đã thử thay đổi thứ tự các bản vá, nhưng tôi không biết cách thực hiện. Hãy cho tôi biết nếu bạn tìm ra cách thực hiện nhé.
## 3.8 Plot the reassigned raster with customized colours# Outline the colours for every class - discover you must comply with the identical order because the values and should be numerically rising 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 = ('Pure Forest', 'Different Pure Vegetation', 'Pasture', 'Agriculture', 'Others', 'Water', 'No knowledge') # 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 top of the listing 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(shade=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 beneath the plot
loc = 'higher middle', # Place the legend within the higher middle
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) # Reserve it as a PDF on the figures folder
plt.present()
4. Các hàm chuẩn hóa
Bây giờ chúng ta đã xây dựng trực giác về cách tải xuống, lưu, dọn dẹp và vẽ các raster LULC của MapBiomas. Đã đến lúc khái quát hóa quá trình.
Trong phần này, chúng tôi sẽ định nghĩa các hàm để tự động hóa các bước này cho bất kỳ hình dạng và năm nào. Sau đó, chúng tôi sẽ thực hiện các hàm này trong một vòng lặp để phân tích một thành phố cụ thể — Porto Acre, AC — từ năm 1985 đến năm 2022. Cuối cùng, chúng tôi sẽ tạo một video minh họa sự phát triển của LULC trong khu vực đó trong khoảng thời gian đã chỉ định.
Đầu tiên, lưu một hộp giới hạn (bbox) xung quanh Vùng quan tâm (ROI). Bạn chỉ cần nhập hình học mong muốn và chỉ định năm. Chức năng này sẽ lưu raster bbox xung quanh ROI vào Google Drive của bạn.
## 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 wish to save the file
fileNamePrefix=f"{geom_name}_{yr}", # File identify
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 venture 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)
Thứ hai, cắt ảnh raster để chỉ bao gồm các điểm ảnh trong hình học và lưu nó dưới dạng ảnh raster mới.
Tôi đã chọn lưu nó trên Google Drive, nhưng bạn có thể thay đổi `reassigned_raster_path` để lưu nó ở bất kỳ nơi nào khác. Nếu bạn thay đổi nó, hãy đảm bảo cập nhật phần còn lại của mã cho phù hợp.
Ngoài ra, bạn có thể chỉ định lại các pixel khi cần bằng cách sửa đổi mapbiomas_categories dict. Chữ số bên trái biểu thị giá trị pixel ban đầu và chữ số bên phải biểu thị giá trị được gán lại (mới).
## 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 may change this to no matter categorization you need, however simply bear 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 exterior 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({
"peak": 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 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 = Exchange, False = Do not substitute)
modified_raster(masks) = new_value # Exchange 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)
Bây giờ chúng ta thấy tần số của từng điểm ảnh trong raster được cắt và chỉ định lại.
## 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)
rely = np.sum(raster_data == worth) # Depend 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}, Depend: {rely}, Share: {share}")
teste3 = pixel_freq_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive)
Cuối cùng, chúng ta vẽ nó trên bản đồ. Bạn có thể thay đổi các đối số bên dưới để điều chỉnh các đặc điểm như màu sắc, nhãn, vị trí chú giải, kích thước phông chữ, v.v. Ngoài ra, còn có tùy chọn để chọn định dạng mà bạn muốn lưu dữ liệu (thường là PDF hoặc PNG). PDF nặng hơn và giữ nguyên độ phân giải, trong khi PNG nhẹ hơn nhưng có độ phân giải thấp hơn.
## 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/[email protected]/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 must comply with 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 = ('Pure Forest', 'Different Pure Vegetation', 'Pasture', 'Agriculture', 'Others', 'Water', 'No knowledge') # Set your labels
cmap = colours.ListedColormap(colors_list) # Create a colormap (cmap) with the colours
bounds = values + (256) # Add a price to the top of the listing 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(shade=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 beneath the plot
loc = 'higher middle', # Place the legend within the higher middle
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) # Reserve it as a .pdf or .png on the figures folder of your venture
plt.present()
teste4 = plot_mapbiomas_lulc_raster('rio_de_janeiro', 2022, folder_in_google_drive, 'png')
Cuối cùng, đây là một ví dụ về cách sử dụng các hàm và tạo vòng lặp để có được sự phát triển của LULC cho Porto Acre (AC) từ năm 1985. Đó là một thành phố khác trong khu vực AMACRO có tỷ lệ phá rừng tăng cao.
## 4.5 Do it in only one perform - 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 identify
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 another shapefile right here
brazilian_municipalities = brazilian_municipalities.clean_names()
brazilian_municipalities.crs = 'EPSG:4326' # Set the CRS to WGS84 (this venture 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 listing 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 the town
geom_name = unidecode.unidecode(metropolis) # Take away latin-1 characters from the town identify - GEE doesn`t enable 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 - Be sure to 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.accumulate()
Chúng tôi sẽ kết thúc hướng dẫn bằng cách tạo một video ngắn cho thấy sự phát triển của nạn phá rừng ở thành phố trong bốn thập kỷ qua. Lưu ý rằng bạn có thể mở rộng phân tích sang nhiều thành phố và chọn những năm cụ thể để phân tích. Hãy thoải mái tùy chỉnh thuật toán khi cần.
## 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.be a part of(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 - ensure you solely have the specified pictures 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)
[ad_2]
Source link