LoginSignup
1
0
@ishimotojuri(Ishimoto Juri)

Google Earth Engineで建物の割合を計算してみる

Posted at

1. はじめに

事業が対象地域の経済活動に貢献しているか検証する際に、建物の割合が一つの指標になりえます。

今回は、Google Earth Engine（GEE）上で建物の割合を計算する方法について紹介したいと思います。

2. データ

今回は、すでに土地利用が分類されている衛星データを使います。
自分で機械学習を使って判別する方法もありますが、より簡便な方法を採用したいと思います。

土地利用データには、Dynamic World V1を使用します。

Dynamic World は、9つ​​の土地利用の書類とその確率、およびラベル情報を含む、解像度10mのデータセットです。GEE上で、2015年6月27日から2024年7月1日までの画像が入手可能です（2024年7月2日時点）。

3. 分析

以下の手順を踏んで建物の割合を計算したいと思います。

  1. 関心領域を設定
  2. データセットを読み込む
  3. 建物のピクセルを抽出
  4. 建物のピクセル数を計算
  5. 結果を出力

3.1 関心領域を設定

関心領域（Area of Interest、AOI）のシェープファイルを用意します。
作成方法はこちらをご覧ください。

AOIのシェープファイルをGoogle Eart Engine（GEE）にインポートします。
インポートの方法はこちらをご参照ください。

今回は、道路をAOIとします。
インポートした道路データを活用して、道路から半径3kmのバッファーを作成します。
分析には、このバッファーデータを使用します。

var target_road =  ee.FeatureCollection('users/ishimotojuri/target_road_4326');
//ご自身のIDを使用してください
var geometryBuffer = target_road.geometry().buffer({'distance':3000}); 
//バッファーの範囲を変えたい場合は、ここを変えます

// AOIをマップに表示させます
Map.addLayer(geometryBuffer, {color: 'red'}, '[red]target_road_buffer');
Map.centerObject(geometryBuffer, 10);

3.2 データセットを読み込む

var collection = ee.ImageCollection("GOOGLE/DYNAMICWORLD/V1") 
            .filterDate("2024-01-01", "2024-06-30") //期間を変えたい場合はここを変えます
            .filterBounds(geometryBuffer)
            .select('label')
            .reduce(ee.Reducer.mode()); 
            //設定した期間の中で最頻値を取得します

3.3 建物のピクセルを抽出

var builtArea = collection.eq(6);
// 0	#419bdf	water
// 1	#397d49	trees
// 2	#88b053	grass
// 3	#7a87c6	flooded_vegetation
// 4	#e49635	crops
// 5	#dfc35a	shrub_and_scrub
// 6	#c4281b	built
// 7	#a59b8f	bare
// 8	#b39fe1	snow_and_ice
// 建物以外を抽出したい場合は、eqの後の数字を変えます

3.4 建物のピクセル数を計算

まず建物のピクセル数を計算します。

var builtArea = builtArea.rename(['built_area']);

// Mask 0 pixel values and count remaining pixels
var builtAreaMasked = builtArea.selfMask();

var statsMasked = builtAreaMasked.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometryBuffer,
    scale: 10,
    maxPixels: 1e10
    });
// "reduceRegion"が、各ピクセルを集計して、任意の領域での統計値
//（ピクセル数、平均値、合計値、最頻値、中央値など）を得るためのコマンドです。

var builtAreaPixels = statsMasked.get('built_area');
print('BuiltArea pixels:', builtAreaPixels);

建物の割合を計算するために、全体のピクセル数も計算します。

// Count totalArea pixels
var composite = collection.rename(['classification']);

var statsTotal = composite.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometryBuffer,
    scale: 10,
    maxPixels: 1e10
    });

var totalPixels = statsTotal.get('classification');
print('Total pixels:', totalPixels);

ここまでで建物の割合は計算できますが、他の土地利用の割合も計算したいという時のために、各土地利用のピクセル数も計算したいと思います。

//Get all pixels by category
var pixelCountStats = composite.reduceRegion({
    reducer: ee.Reducer.frequencyHistogram().unweighted(),
    geometry: geometryBuffer,
    scale: 10,
    maxPixels: 1e10
    });


var pixelCounts = ee.Dictionary(pixelCountStats.get('classification'));

// Format the results to make it readable
// var classLabels = ee.List([
//     'water', 'trees', 'grass', 'flooded_vegetation', 'crops',
//     'shrub_and_scrub', 'built', 'bare', 'snow_and_ice'
//     ]);


// There is no snow_and_ice in 2018, so we need to exclude it from the list.
var classLabels = ee.List([
    'water', 'trees', 'grass', 'flooded_vegetation', 'crops',
    'shrub_and_scrub', 'built', 'bare'
    ]);
// 年によってはすべての土地利用のピクセルがないこともあります。
// その場合、リストに含めてしまっているとエラーがでてしまうので、適宜リストから除外します。
// 'snow_and_ice'がない場合が多いです。上のリストでは除外しています。

// Rename keys with class names
var pixelCountsFormatted = pixelCounts.rename(pixelCounts.keys(), classLabels);

print('Total pixels by category:', pixelCountsFormatted);

3.5 結果を出力

上記で計算した各土地利用のピクセル数をエクセル形式でダウンロードします。

// Create a Feature Collection
var exportFc = ee.FeatureCollection(ee.Feature(null, pixelCountsFormatted));

// Export the results as a CSV file
Export.table.toDrive({
  collection: exportFc,
  description: 'pixel_counts_export_1km', //change here 
  folder: 'Vietnam',
  fileNamePrefix: 'pixel_counts_1km', //change here 
  fileFormat: 'CSV',
});

4. まとめ

以下、一連のコードです。

////////////////////////////////////////////////////////////////////////////////
// 1. Set an AOI
////////////////////////////////////////////////////////////////////////////////
var target_road =  ee.FeatureCollection('users/ishimotojuri/target_road_4326');

var geometryBuffer = target_road.geometry().buffer({'distance':3000}); //change here

Map.addLayer(geometryBuffer, {color: 'red'}, '[red]target_road_buffer');
Map.centerObject(geometryBuffer, 10);

// var londuc = ee.FeatureCollection('users/ishimotojuri/londuc');
// Map.addLayer(londuc, {color: 'FF0000'}, 'londuc');

// var urban_area = ee.FeatureCollection('users/ishimotojuri/urban');
// Map.addLayer(urban_area, {color:'999999'}, 'urban_area');

// var catlai = ee.FeatureCollection('users/ishimotojuri/CatLai');
// Map.addLayer(catlai, {color: 'orange'}, 'catlai');


////////////////////////////////////////////////////////////////////////////////
// 2. Import the dataset
///////////////////////////////////////////////////////////////////////////////
var collection = ee.ImageCollection("GOOGLE/DYNAMICWORLD/V1") 
            .filterDate("2024-01-01", "2024-06-30") //Change here
            .filterBounds(geometryBuffer)
            .select('label')
            .reduce(ee.Reducer.mode()); 
            //最頻値

////////////////////////////////////////////////////////////////////////////////
// 3. Extract builtAreas
////////////////////////////////////////////////////////////////////////////////

var builtArea = collection.eq(6);
// 0	#419bdf	water
// 1	#397d49	trees
// 2	#88b053	grass
// 3	#7a87c6	flooded_vegetation
// 4	#e49635	crops
// 5	#dfc35a	shrub_and_scrub
// 6	#c4281b	built
// 7	#a59b8f	bare
// 8	#b39fe1	snow_and_ice

////////////////////////////////////////////////////////////////////////////////
// 4. Count builtArea pixels
////////////////////////////////////////////////////////////////////////////////
var builtArea = builtArea.rename(['built_area']);

// Mask 0 pixel values and count remaining pixels
var builtAreaMasked = builtArea.selfMask();

var statsMasked = builtAreaMasked.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometryBuffer,
    scale: 10,
    maxPixels: 1e10
    });

var builtAreaPixels = statsMasked.get('built_area');
print('BuiltArea pixels:', builtAreaPixels);

// Count totalArea pixels
var composite = collection.rename(['classification']);

var statsTotal = composite.reduceRegion({
    reducer: ee.Reducer.count(),
    geometry: geometryBuffer,
    scale: 10,
    maxPixels: 1e10
    });

var totalPixels = statsTotal.get('classification');
print('Total pixels:', totalPixels);


//Get all pixels by category
var pixelCountStats = composite.reduceRegion({
    reducer: ee.Reducer.frequencyHistogram().unweighted(),
    geometry: geometryBuffer,
    scale: 10,
    maxPixels: 1e10
    });

var pixelCounts = ee.Dictionary(pixelCountStats.get('classification'));

// Format the results to make it readable
// var classLabels = ee.List([
//     'water', 'trees', 'grass', 'flooded_vegetation', 'crops',
//     'shrub_and_scrub', 'built', 'bare', 'snow_and_ice'
//     ]);


// There is no snow_and_ice in 2018, so we need to exclude it from the list.
var classLabels = ee.List([
    'water', 'trees', 'grass', 'flooded_vegetation', 'crops',
    'shrub_and_scrub', 'built', 'bare',
    ]);


// Rename keys with class names
var pixelCountsFormatted = pixelCounts.rename(pixelCounts.keys(), classLabels);

print('Total pixels by category:', pixelCountsFormatted);

////////////////////////////////////////////////////////////////////////////////
//5. Export the result
////////////////////////////////////////////////////////////////////////////////

// Create a Feature Collection
var exportFc = ee.FeatureCollection(ee.Feature(null, pixelCountsFormatted));

// Export the results as a CSV file
Export.table.toDrive({
  collection: exportFc,
  description: 'pixel_counts_export_1km', //change here 
  folder: 'Vietnam',
  fileNamePrefix: 'pixel_counts_1km', //change here 
  fileFormat: 'CSV',
});

5. おわりに

今回は建物の割合の計算方法をご紹介しました。

期間を変えれば、複数年の建物の割合を計算できるため、時系列分析もできます。

ご関心の領域で試してみてください！

ソースコード

  • 一連のコードはこちらからもご覧いただけます。
  • 道路のシェープファイルはご自身でインポートし、Table IDを再設定する必要があるので、ご注意ください。
  • コードの解釈や実行に関する質問があれば、お気軽にコメントしてください！
1
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
1
0