Первым шагом в классификации землепользования является получение данных. Это можно сделать с помощью спутниковых снимков, аэрофотоснимков или наземных съемок. После того, как данные получены, их необходимо предварительно обработать, чтобы удалить любой шум или артефакты, которые могут помешать процессу классификации. Это может включать удаление облаков или других атмосферных эффектов и корректировку различий в освещении. В этом уроке я использовал пример изображения Sentinel-2 над Токио, Япония, с полосами [синий, зеленый, красный и NIR] с минимальной облачностью. Не стесняйтесь использовать другие доступные источники данных.
Спутниковое изображение. Вы можете загрузить любое спутниковое изображение из USGS или GEE.
Ярлыки (образцы). Вы можете собрать информацию либо из полевого опроса, либо вручную с помощью наложения изображения. Для ручного создания используйте любое программное обеспечение ГИС, создайте шейп-файл (точку) и добавьте точки с соответствующими метками в качестве атрибутов.
Я загрузил входное изображение и метки в этот репозиторий блога GitHub (ссылка в конце). Смело используйте их сами!
После того, как вы завершили загрузку спутникового снимка и собрали метки (образцы), мы можем перейти к рабочему процессу кодирования.
2. Извлечение признаков:
Во-первых, мы импортируем необходимые библиотеки Python, предоставим входные местоположения для нашего изображения и меток, а также список имен каналов входного изображения.
import rasterio as rio from rasterio.plot import show import geopandas as gpd import fiona import matplotlib.pyplot as plt import seaborn as sns import numpy as np import pandas as pd import sklearn from sklearn.metrics import classification_report, accuracy_score # variables # Note: labels should be always last column with name "labels" # Note: Make sure input labels shapefile and input raster have same CRS, otherwise code will not run # input files raster_loc = 'materials/rasters/s2image.tif' points_loc = 'materials/shapefiles/samples.shp' temp_point_loc = 'materials/temp/temp_y_points.shp' # land cover names (for post visualization) lulc_name = ['Water', 'Dense Veg', 'Veg', 'Impervious']
Необязательный шаг: вы можете визуализировать доступные каналы спутниковых изображений вместе с композициями естественного цвета и искусственного цвета (коды приведены в блокноте Jupyter), используя подграфики matplotlib. В результате полосы будут выглядеть так:
Затем прочитайте входное спутниковое изображение как растерический массив и пометьте шейп-файл как кадр данных geopandas. Для извлечения признаков сначала присвойте уникальный идентификационный номер и сохраните его как новый временный файл.
# reading bands from input
with rio.open(raster_loc) as img:
bands = (img.read()).shape[0]
print('Bands of input image: ', bands)
# using ilteration to automatically create a bands list
features = []
for i in range(bands):
features.append('band'+str(i+1))
print('Bands names: ', features)
f_len = len(features)
points = gpd.read_file(points_loc)
# adding a new column 'id' with range of points
points = points.assign(id=range(len(points)))
# saving nenw point file with 'id'
points.to_file(temp_point_loc)
# converting gdf to pd df and removing geometry
points_df = pd.DataFrame(points.drop(columns='geometry'))
Затем создайте пустой кадр данных серии pandas и повторите каждую метку, используя Fiona с циклом python for. Для каждой итерации соответствующее значение каждой полосы (и других функций, таких как NDVI) сохраняется для каждого местоположения метки. После итерации измените форму данных в формате списка.
# ilterating over multiband raster
sampled = pd.Series()
#inputShape= temp_point_loc
# Read input shapefile with fiona and iterate over each feature
with fiona.open(temp_point_loc) as shp:
for feature in shp:
siteID = feature['properties']['id']
coords = feature['geometry']['coordinates']
# Read pixel value at the given coordinates using Rasterio
# NB: `sample()` returns an iterable of ndarrays.
with rio.open(raster_loc) as stack_src:
value = [v for v in stack_src.sample([coords])]
# Update the pandas serie accordingly
sampled.loc[siteID] = value
# reshaping sampled values
df1 = pd.DataFrame(sampled.values.tolist(), index=sampled.index)
df1['id'] = df1.index
df1 = pd.DataFrame(df1[0].values.tolist(),
columns=features)
df1['id'] = df1.index
data = pd.merge(df1, points_df, on ='id')
print('Sampled Data: \n',data)
После выборки признаков разделите данные на X (независимая переменная — признаки) и Y (зависимая переменная — метки). После этого используйте функцию Scikit-learn train_test_split, чтобы разделить данные в соотношении 70/30 (70% для обучения и 30% для тестирования).
x = data.iloc[:,0:f_len]
X = x.values
y = data.iloc[:,-1]
Y = y.values
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, Y, test_size=0.30, stratify = Y)
print(f'X_train Shape: {X_train.shape}\nX_test Shape: {X_test.shape}\ny_train Shape: {y_train.shape}\ny_test Shape:{y_test.shape}')
3. Обучение модели и оценка точности
После подготовки данных для обучения и тестирования очень просто запустить модели с параметрами по умолчанию (тема Гипараметрическая оптимизация может быть освещена в другом блоге).
SVM
Во-первых, мы начнем с модели машины опорных векторов (SVM). Сначала импортируйте модель (например, из Scikit-learn). Затем определите его параметры (в случае SVM мы вводим ядро как rbf). После инициализации модели используйте функцию fit для обучения модели с использованием данных x и y. После этого используйте тестовые данные x, чтобы проверить прогноз выборки, и используйте Scikit-learn accuracy_score, чтобы проверить точность модели для этого конкретного прогноза выборки. В нашем случае SVM правильно классифицировал все пиксели, поэтому точность 100%.
from sklearn.svm import SVC
from sklearn.metrics import confusion_matrix
cName = 'SVM'
clf = SVC(kernel='rbf')
clf.fit(X_train, y_train)
clf_pred = clf.predict(X_test)
print(f"Accuracy {cName}: {accuracy_score(y_test, clf_pred)*100}")
print(classification_report(y_test, clf_pred))
Чтобы дополнительно проверить производительность модели, мы можем оценить матрицу путаницы и визуализировать ее с помощью функции Seaborn ‘‘sns.heatmap’’.
# Confusion Matrix
cm = confusion_matrix(y_test, clf_pred)
print('Confusion Matrix RF: \n',cm)
cm_percent = cm/np.sum(cm)
plt.figure(figsize=(7, 7), facecolor='w', edgecolor='k')
sns.set(font_scale=1.5)
sns.heatmap(cm_percent,
xticklabels=lulc_name,
yticklabels=lulc_name,
cmap="YlGn",
annot=True,
fmt='.2%',
cbar=False,
linewidths=2,
linecolor='black')
plt.title(cName)
plt.xlabel('Predicted')
plt.ylabel('Actual')
#plt.savefig(f'../figs/confusion_matrix_{cName}.png', dpi=300, bbox_inches='tight')
Результирующая матрица путаницы выглядит следующим образом. Чтобы узнать о дальнейших изменениях, см. этот блог.
Теперь мы можем выполнить те же действия для разных моделей. Список доступных моделей классификации в Scikit-learn смотрите в этом списке. Вы также можете использовать другие модели (вне Scikit-learn), такие как XGBoost и LightGBM. Для дальнейшего тестирования мы запускаем модели Random Forest (RF) и Decision Trees (DT), матрица путаницы которых приведена ниже (коды доступны в основном блокноте Jupyter).
RF имеет точность 94%, тогда как DT имеет точность 89%. Так что в нашем случае и SVM, и RF показали хорошие результаты. Также обратите внимание, что максимальная путаница наблюдается между растительностью и непроницаемым классом LULC.
Есть и другие способы оценки производительности модели классификации, такие как кривая ROC/AUC, важность переменной модели, точность отзыва, f1-баллы и т. д. Я расскажу о них в своих будущих блогах.
4. Прогнозирование и экспорт данных
Наиболее важным шагом является прогнозирование полных данных с использованием нашей обученной модели и экспорт численно прогнозируемых результатов в собственный формат Geotiff. Хотя доступно много вариантов (GDAL или использование GeoPandas/Geocube), они гораздо менее эффективны при масштабировании этого процесса на более крупные географические регионы (в настоящее время я занимаюсь анализом в национальном масштабе, и эти методы являются узкими местами). Таким образом, лучший способ использования rasterio дается следующим образом.
Сначала прочитайте исходное входное изображение и получите различные свойства метаданных (высота, ширина, CRS). После этого измените форму изображения так же, как вы предоставили входные данные для обучения. Обратите внимание: если вы ранее добавили дополнительные функции, такие как NDVI или Elevation, вам может потребоваться сначала скомпоновать их с входным изображением, а затем выполнить этот процесс. Используйте для этого код ниже:
%%time
cName = 'SVM'
exp_name = f'materials/results/lulc_{cName}.tif'
img = rio.open(raster_loc)
img_arr = img.read()
bands = img_arr.shape[0]
print(f'Height: {img_arr.shape[1]}\nWidth: {img_arr.shape[2]}\nBands: {img_arr.shape[0]}\n')
img_n = np.moveaxis(img_arr, 0, -1)
img_n = img_n.reshape(-1, f_len)
print('reshaped full data shape for prediction: ',img_n.shape)
metadata = img.met
height = metadata.get('height')
width = metadata.get('width')
crs = metadata.get('crs')
transform = metadata.get('transform')
Затем используйте существующую модель с командой прогнозирования, чтобы предсказать полные данные, которые вы только что изменили. В зависимости от ваших входных данных, это процесс, который займет большую часть времени (и использование вашей оперативной памяти и скорости процессора). После предсказания снова измените его форму, чтобы ввести (исходная высота, ширина изображения).
pred_full = clf.predict(img_n)
print('Prediction Done, now exporting raster \n')
img_reshape = pred_full.reshape(height, width)
Наконец, используйте команду rasterio export raster для экспорта предсказанных/измененных результатов.
out_raster = rio.open(exp_name,
'w',
driver='GTiff',
height=height,
width=width,
count=1, # output band number
dtype='uint8', #output data type
crs=crs,
transform = transform,
nodata = 255 #nodata
)
out_raster.write(img_reshape, 1)
out_raster.close()
После экспорта вы можете открыть его с помощью Python или программного обеспечения Географической информационной системы (ГИС) (ArcGIS Pro или QGIS). Вывод будет выглядеть так (конечно, после того, как вы установите цвета, имена и ключ карты и т. д.)
Вы также можете наложить свой LULC на входное изображение и посмотреть, соответствуют ли классы правильным меткам/местоположениям.
Заключение:
В этом блоге мы рассмотрели, как классифицировать LULC на основе ML с помощью Python. Процесс включает в себя сбор и предварительную обработку данных, извлечение признаков, обучение модели классификации и классификацию новых данных. Python имеет несколько библиотек, которые можно использовать для каждого из этих шагов, что делает его мощным инструментом для классификации землепользования. При наличии правильных инструментов и методов классификация землепользования может помочь нам лучше понять наши природные ресурсы и управлять ими. Коды, связанные с этим блогом, находятся в свободном доступе на моем GitHub (https://github.com/waleedgeo/lulc_py).
Если этот блог был вам полезен, поделитесь им и подпишитесь на меня, чтобы получать больше такого контента. Кроме того, вы можете использовать мои опубликованные исследования для справки по LULC на основе ML.
Мой LinkedIn: https://www.linkedin.com/in/waleedgeo/
Веб-страница моего портфолио: www.waleedgeo.com
— — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — — —
Благодарность
Этот блог был бы невозможен без вклада с открытым исходным кодом:
Люди:
Доктор. Саджад Мухаммад (мой научный руководитель, который дал мне эту идею для работы)
Проф. Qiusheng Wu (Бесчисленный вклад в сообщество специалистов по геопространственным данным с открытым исходным кодом). Проверьте его пакеты Geemap и Leafmap.
Бонни П (автор книги Python для анализа геопространственных данных)
Блоги:
Книги:
Прикладное прогнозное моделирование Макса Куна
Руководство для студентов по Python для физического моделирования, Джесси М. Киндер
