การแยกประเภทแบบไม่มีการควบคุมด้วยชุดข้อมูลการฝังดาวเทียม

แก้ไขใน GitHub
รายงานปัญหา
ประวัติหน้าเว็บ
ผู้เขียน spatialthoughts
บทแนะนำนี้เป็นส่วนหนึ่งของชุดบทแนะนำเกี่ยวกับชุดข้อมูลการฝังดาวเทียม โปรดดู ข้อมูลเบื้องต้น การจัดประเภทภายใต้การกำกับดูแล การถดถอย และการค้นหาความคล้ายคลึงกันด้วย

ในบทแนะนำก่อนหน้า (ข้อมูลเบื้องต้น) เราได้เห็นว่าการฝังดาวเทียมบันทึกวิถีประจำปีของการสังเกตการณ์ดาวเทียมและตัวแปรสภาพอากาศได้อย่างไร ซึ่งทำให้ชุดข้อมูลพร้อมใช้งานสำหรับการทำแผนที่พืชผลโดยไม่ต้องสร้างแบบจำลองฟีโนโลยีของพืชผล การแมปประเภทพืชเป็นงานที่ท้าทายซึ่งโดยปกติแล้วต้องมีการสร้างแบบจำลองปรากฏการณ์ทางสรีรวิทยาของพืชและการเก็บตัวอย่างในแปลงสำหรับพืชทั้งหมดที่ปลูกในภูมิภาค

ในบทแนะนำนี้ เราจะใช้แนวทางการแยกประเภทแบบไม่มีการกำกับดูแลในการทำแผนที่พืชผล ซึ่งจะช่วยให้เราทำงานที่ซับซ้อนนี้ได้โดยไม่ต้องอาศัยป้ายกำกับฟิลด์ วิธีนี้ใช้ประโยชน์จากความรู้ในท้องถิ่นของภูมิภาคร่วมกับสถิติพืชผลรวม ซึ่งพร้อมให้บริการในหลายส่วนของโลก

เลือกภูมิภาค

ในบทแนะนำนี้ เราจะสร้างแผนที่ประเภทพืชสำหรับเคาน์ตีเซร์โร กอร์โด รัฐไอโอวา เขตนี้อยู่ในแถบปลูกข้าวโพดของสหรัฐอเมริกา ซึ่งมีพืชหลัก 2 ชนิด ได้แก่ ข้าวโพดและถั่วเหลือง ความรู้ในท้องถิ่นนี้มีความสำคัญและจะช่วยให้เราตัดสินใจเกี่ยวกับพารามิเตอร์ที่สำคัญสำหรับโมเดลของเราได้

มาเริ่มจากการขอขอบเขตสำหรับเขตที่เลือกกัน

// Select the region
// Cerro Gordo County, Iowa
var counties = ee.FeatureCollection('TIGER/2018/Counties');

var selected = counties
  .filter(ee.Filter.eq('GEOID', '19033'));
var geometry = selected.geometry();
Map.centerObject(geometry, 12);
Map.addLayer(geometry, {color: 'red'}, 'Selected Region', false);


รูปภาพ: ภูมิภาคที่เลือก

เตรียมชุดข้อมูลการฝังดาวเทียม

จากนั้นเราจะโหลดชุดข้อมูลการฝังดาวเทียม กรองรูปภาพสำหรับปีที่เลือก และสร้างภาพโมเสก

var embeddings = ee.ImageCollection('GOOGLE/SATELLITE_EMBEDDING/V1/ANNUAL');

var year = 2022;
var startDate = ee.Date.fromYMD(year, 1, 1);
var endDate = startDate.advance(1, 'year');

var filteredembeddings = embeddings
  .filter(ee.Filter.date(startDate, endDate))
  .filter(ee.Filter.bounds(geometry));

var embeddingsImage = filteredembeddings.mosaic();

สร้างมาสก์พื้นที่เพาะปลูก

สำหรับการประมาณ เราต้องยกเว้นพื้นที่ที่ไม่ใช่พื้นที่เพาะปลูก มีชุดข้อมูลระดับโลกและระดับภูมิภาคมากมายที่ใช้สร้างมาสก์ครอบตัดได้ ESA WorldCover หรือ GFSAD Global Cropland Extent Product เป็นตัวเลือกที่ดีสำหรับชุดข้อมูลพื้นที่เพาะปลูกทั่วโลก นอกจากนี้ เรายังได้เพิ่มผลิตภัณฑ์ ESA WorldCereal Active Cropland ซึ่งมีการแมปพื้นที่เพาะปลูกที่ใช้งานอยู่ตามฤดูกาล เนื่องจากภูมิภาคของเราอยู่ในสหรัฐอเมริกา เราจึงใช้ชุดข้อมูลระดับภูมิภาคที่แม่นยำยิ่งขึ้นอย่างเลเยอร์ข้อมูลพื้นที่เพาะปลูกของ NASS จาก USDA (CDL) เพื่อรับมาสก์พืชผลได้

// Use Cropland Data Layers (CDL) to obtain cultivated cropland
var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var croplandMask = cdl.select('cultivated').eq(2).rename('cropmask');

// Visualize the crop mask
var croplandMaskVis = {min: 0, max: 1, palette: ['white', 'green']};
Map.addLayer(croplandMask.clip(geometry), croplandMaskVis, 'Crop Mask');


รูปภาพ: ภูมิภาคที่เลือกพร้อมมาสก์พื้นที่เพาะปลูก

แยกตัวอย่างการฝึก

เราใช้มาสก์พื้นที่เพาะปลูกกับโมเสกการฝัง ตอนนี้เราเหลือเพียงพิกเซลทั้งหมดที่แสดงถึงพื้นที่เพาะปลูกในเขต

// Mask all non-cropland pixels
var clusterImage = embeddingsImage.updateMask(croplandMask);

เราต้องนำรูปภาพการฝังดาวเทียมและสุ่มตัวอย่างเพื่อฝึกโมเดลการจัดกลุ่ม เนื่องจากภูมิภาคที่เราสนใจมีพิกเซลที่มาสก์ไว้จำนวนมาก การสุ่มตัวอย่างแบบสุ่มอย่างง่ายอาจทำให้ได้ตัวอย่างที่มีค่าเป็น Null เราใช้การสุ่มตัวอย่างแบบแบ่งชั้นเพื่อให้แน่ใจว่าเราสามารถดึงตัวอย่างที่ไม่ใช่ค่าว่างตามจำนวนที่ต้องการได้ เพื่อให้ได้ตัวอย่างตามจำนวนที่ต้องการในพื้นที่ที่ไม่ได้มาสก์

// Stratified random sampling
var training = clusterImage.addBands(croplandMask).stratifiedSample({
  numPoints: 1000,
  classBand: 'cropmask',
  region: geometry,
  scale: 10,
  tileScale: 16,
  seed: 100,
  dropNulls: true,
  geometries: true
});

ส่งออกตัวอย่างไปยังชิ้นงาน (ไม่บังคับ)

การแยกตัวอย่างเป็นการดำเนินการที่ใช้การคำนวณสูง และแนวทางปฏิบัติที่ดีคือการส่งออกตัวอย่างการฝึกที่แยกออกมาเป็นชิ้นงาน แล้วใช้ชิ้นงานที่ส่งออกในขั้นตอนต่อๆ ไป ซึ่งจะช่วยแก้ไขข้อผิดพลาดการคำนวณหมดเวลาหรือหน่วยความจำของผู้ใช้เกินเมื่อทำงานกับภูมิภาคขนาดใหญ่

เริ่มงานส่งออกและรอให้เสร็จสิ้นก่อนดำเนินการต่อ

// Replace this with your asset folder
// The folder must exist before exporting
var exportFolder = 'projects/spatialthoughts/assets/satellite_embedding/';

var samplesExportFc = 'cluster_training_samples';
var samplesExportFcPath = exportFolder + samplesExportFc;

Export.table.toAsset({
  collection: training,
  description: 'Cluster_Training_Samples',
  assetId: samplesExportFcPath
});

เมื่องานส่งออกเสร็จสมบูรณ์แล้ว เราจะอ่านตัวอย่างที่แยกออกมากลับเข้าไปในโค้ดเป็นคอลเล็กชันฟีเจอร์ได้

// Use the exported asset
var training = ee.FeatureCollection(samplesExportFcPath);

แสดงตัวอย่างเป็นภาพ

ไม่ว่าคุณจะเรียกใช้ตัวอย่างแบบอินเทอร์แอกทีฟหรือส่งออกไปยังคอลเล็กชันฟีเจอร์ ตอนนี้คุณจะมีตัวแปรการฝึกพร้อมจุดตัวอย่าง มาลองพิมพ์ตัวอย่างแรกเพื่อตรวจสอบและเพิ่มจุดฝึกไปยัง Map กัน

print('Extracted sample', training.first());
Map.addLayer(training, {color: 'blue'}, 'Extracted Samples', false);


รูปภาพ: ตัวอย่างแบบสุ่มที่ดึงออกมาสำหรับการจัดกลุ่ม

ทำการจัดกลุ่มแบบไม่มีการควบคุมดูแล

ตอนนี้เราสามารถฝึกเครื่องจัดกลุ่มและจัดกลุ่มเวกเตอร์การฝัง 64 มิติเป็นคลัสเตอร์ที่แตกต่างกันตามจำนวนที่เลือกได้แล้ว จากความรู้ในท้องถิ่น เราพบว่ามีพืชผล 2 ประเภทหลักที่ครอบคลุมพื้นที่ส่วนใหญ่ และมีพืชผลอื่นๆ อีกหลายชนิดที่ครอบคลุมพื้นที่ส่วนที่เหลือ เราสามารถทำการจัดกลุ่มแบบไม่มีการกำกับดูแลในข้อมูลฝังดาวเทียมเพื่อรับกลุ่มพิกเซลที่มีวิถีและรูปแบบชั่วคราวที่คล้ายกัน พิกเซลที่มีลักษณะทางสเปกตรัมและเชิงพื้นที่คล้ายกัน รวมถึงฟีโนโลยีที่คล้ายกันจะได้รับการจัดกลุ่มไว้ในคลัสเตอร์เดียวกัน

ee.Clusterer.wekaCascadeKMeans() ช่วยให้เราสามารถระบุจำนวนคลัสเตอร์ขั้นต่ำและสูงสุด รวมถึงค้นหาจำนวนคลัสเตอร์ที่เหมาะสมที่สุดตามข้อมูลการฝึก ความรู้ในพื้นที่ของเราจึงมีประโยชน์ในการกำหนดจำนวนคลัสเตอร์ขั้นต่ำและสูงสุด เนื่องจากเราคาดว่าจะมีคลัสเตอร์ที่แตกต่างกัน 2-3 ประเภท ได้แก่ ข้าวโพด ถั่วเหลือง และพืชอื่นๆ อีกหลายชนิด เราจึงใช้ 4 เป็นจำนวนคลัสเตอร์ขั้นต่ำ และ 5 เป็นจำนวนคลัสเตอร์สูงสุด คุณอาจต้องทดลองกับตัวเลขเหล่านี้เพื่อดูว่าอะไรเหมาะกับภูมิภาคของคุณมากที่สุด

// Cluster the Satellite Embedding Image
var minClusters = 4;
var maxClusters = 5;

var clusterer = ee.Clusterer.wekaCascadeKMeans({
  minClusters: minClusters, maxClusters: maxClusters}).train({
  features: training,
  inputProperties: clusterImage.bandNames()
});

var clustered = clusterImage.cluster(clusterer);
Map.addLayer(clustered.randomVisualizer().clip(geometry), {}, 'Clusters');


รูปภาพ: คลัสเตอร์ที่ได้จากการจัดประเภทแบบไม่มีการกำกับดูแล

กำหนดป้ายกำกับให้กับคลัสเตอร์

จากการตรวจสอบด้วยสายตา พบว่ากลุ่มที่ได้ในขั้นตอนก่อนหน้าตรงกับขอบเขตฟาร์มที่เห็นในรูปภาพความละเอียดสูง เราทราบจากความรู้ในท้องถิ่นว่ากลุ่มที่ใหญ่ที่สุด 2 กลุ่มคือข้าวโพดและถั่วเหลือง มาคำนวณพื้นที่ของแต่ละคลัสเตอร์ในรูปภาพกัน

// Calculate Cluster Areas
// 1 Acre = 4046.86 Sq. Meters
var areaImage = ee.Image.pixelArea().divide(4046.86).addBands(clustered);

var areas = areaImage.reduceRegion({
      reducer: ee.Reducer.sum().group({
      groupField: 1,
      groupName: 'cluster',
    }),
    geometry: geometry,
    scale: 10,
    maxPixels: 1e10
    });

var clusterAreas = ee.List(areas.get('groups'));

// Process results to extract the areas and create a FeatureCollection

var clusterAreas = clusterAreas.map(function(item) {
  var areaDict = ee.Dictionary(item);
  var clusterNumber = areaDict.getNumber('cluster').format();
  var area = areaDict.getNumber('sum')
  return ee.Feature(null, {cluster: clusterNumber, area: area})
})

var clusterAreaFc = ee.FeatureCollection(clusterAreas);
print('Cluster Areas', clusterAreaFc);

เราจะเลือก 2 คลัสเตอร์ที่มีพื้นที่สูงสุด

var selectedFc = clusterAreaFc.sort('area', false).limit(2);
print('Top 2 Clusters by Area', selectedFc);

แต่เรายังไม่ทราบว่าคลัสเตอร์ใดเป็นพืชชนิดใด หากมีตัวอย่างฟิลด์ข้าวโพดหรือถั่วเหลือง 2-3 ตัวอย่าง คุณสามารถวางซ้อนบนคลัสเตอร์เพื่อดูป้ายกำกับที่เกี่ยวข้องได้ หากไม่มีตัวอย่างฟิลด์ เราจะใช้ประโยชน์จากสถิติพืชผลรวมได้ ในหลายพื้นที่ทั่วโลก มีการรวบรวมและเผยแพร่สถิติพืชผลรวมเป็นประจำ สำหรับสหรัฐอเมริกา บริการสถิติการเกษตรแห่งชาติ (NASS) มีสถิติพืชผลโดยละเอียดสำหรับแต่ละเทศมณฑลและพืชผลหลักแต่ละชนิด ในปี 2022 เขต Cerro Gordo รัฐไอโอวามีพื้นที่ปลูกข้าวโพด 161,500 เอเคอร์ และพื้นที่ปลูกถั่วเหลือง 110,500 เอเคอร์

เมื่อใช้ข้อมูลนี้ เราจึงทราบว่าในกลุ่ม 2 กลุ่มแรก กลุ่มที่มีพื้นที่ใหญ่ที่สุดน่าจะเป็นข้าวโพด และอีกกลุ่มหนึ่งน่าจะเป็นถั่วเหลือง มากำหนดป้ายกำกับเหล่านี้และเปรียบเทียบพื้นที่ที่คำนวณแล้วกับสถิติที่เผยแพร่กัน

var cornFeature = selectedFc.sort('area', false).first();
var soybeanFeature = selectedFc.sort('area').first();
var cornCluster = cornFeature.get('cluster');
var soybeanCluster = soybeanFeature.get('cluster');

print('Corn Area (Detected)', cornFeature.getNumber('area').round());
print('Corn Area (From Crop Statistics)', 163500);

print('Soybean Area (Detected)', soybeanFeature.getNumber('area').round());
print('Soybean Area (From Crop Statistics)', 110500);

สร้างแผนที่พืช

ตอนนี้เรารู้ป้ายกำกับสำหรับแต่ละคลัสเตอร์แล้ว และสามารถดึงพิกเซลสำหรับแต่ละประเภทการครอบตัดและผสานรวมเพื่อสร้างแผนที่การครอบตัดสุดท้าย

// Select the clusters to create the crop map
var corn = clustered.eq(ee.Number.parse(cornCluster));
var soybean = clustered.eq(ee.Number.parse(soybeanCluster));

var merged = corn.add(soybean.multiply(2));
var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(merged.clip(geometry), cropVis, 'Crop Map (Detected)');

นอกจากนี้ เรายังใช้คอมโพเนนต์ UI เพื่อสร้างและเพิ่มคำอธิบายแทนลงในแผนที่เพื่อช่วยตีความผลลัพธ์ได้ด้วย

// Add a Legend
var legend = ui.Panel({
  layout: ui.Panel.Layout.Flow('horizontal'),
  style: {position: 'bottom-center', padding: '8px 15px'}});

var addItem = function(color, name) {
  var colorBox = ui.Label({
    style: {color: '#ffffff',
      backgroundColor: color,
      padding: '10px',
      margin: '0 4px 4px 0',
    }
  });
  var description = ui.Label({
    value: name,
    style: {
      margin: '0px 10px 0px 2px',
    }
  });
  return ui.Panel({
    widgets: [colorBox, description],
    layout: ui.Panel.Layout.Flow('horizontal')}
)};

var title = ui.Label({
  value: 'Legend',
  style: {fontWeight: 'bold',
    fontSize: '16px',
    margin: '0px 10px 0px 4px'}});

legend.add(title);
legend.add(addItem('#ffd400', 'Corn'));
legend.add(addItem('#267300', 'Soybean'));
legend.add(addItem('#bdbdbd', 'Other Crops'));


รูปภาพ: แผนที่พืชผลที่ตรวจพบซึ่งมีข้าวโพดและถั่วเหลือง

ตรวจสอบผลลัพธ์

เราสามารถรับแผนที่ประเภทพืชผลที่มีชุดข้อมูลการฝังดาวเทียมโดยไม่มีป้ายกำกับฟิลด์ใดๆ โดยใช้เพียงสถิติรวมและความรู้ในท้องถิ่นของภูมิภาค มาเปรียบเทียบผลลัพธ์ของเรากับแผนที่ประเภทพืชอย่างเป็นทางการจากชั้นข้อมูลพื้นที่เพาะปลูก (CDL) ของ NASS จาก USDA กัน

var cdl = ee.ImageCollection('USDA/NASS/CDL')
  .filter(ee.Filter.date(startDate, endDate))
  .first();
var cropLandcover = cdl.select('cropland');
var cropMap = cropLandcover.updateMask(croplandMask).rename('crops');

// Original data has unique values for each crop ranging from 0 to 254
var cropClasses = ee.List.sequence(0, 254);
// We remap all values as following
// Crop     | Source Value | Target Value
// Corn     | 1            | 1
// Soybean  | 5            | 2
// All other| 0-255        | 0
var targetClasses = ee.List.repeat(0, 255).set(1, 1).set(5, 2);
var cropMapReclass = cropMap.remap(cropClasses, targetClasses).rename('crops');

var cropVis = {min: 0, max: 2, palette: ['#bdbdbd', '#ffd400', '#267300']};
Map.addLayer(cropMapReclass.clip(geometry), cropVis, 'Crop Landcover (CDL)');

แม้ว่าผลลัพธ์ของเราจะแตกต่างจากแผนที่อย่างเป็นทางการ แต่คุณจะเห็นว่าเราได้รับผลลัพธ์ที่ดีมากโดยใช้ความพยายามเพียงเล็กน้อย การใช้ขั้นตอนการประมวลผลภายหลังกับผลลัพธ์จะช่วยให้เรานำนอยส์บางส่วนออกและเติมช่องว่างในเอาต์พุตได้

ลองใช้สคริปต์ฉบับเต็มสำหรับบทแนะนำนี้ในโปรแกรมแก้ไขโค้ดของ Earth Engine