GEE8：多个矢量点的NDVI连续数据的获取及分析【CSV数据】

1. 哨兵2号数据：

Sentinel-2是一个宽波段、高分辨率、多光谱成像任务，支持哥白尼土地监测研究，包括监测植被、土壤和水覆盖，以及观察内陆水道和沿海地区。Sentinel-2数据包含13个UINT16波段，代表TOA反射率10000倍。此外，有三个QA波段，其中一个(QA60)是带云掩码信息的位掩码波段。

1.1 Sentinel-2 MSI: MultiSpectral Instrument, Level-1C

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000);
}

// Map the function over one month of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2018-01-01', '2018-01-31')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

Map.setCenter(113, 31, 7);
Map.addLayer(dataset.median(), rgbVis, 'RGB');

天津市图像：

1.2 GEE获取矢量站点连续NDVI值

1.2.1 数据导入

var p1 = ee.Geometry.Point([112.91, 28.24])
var p2 = ee.Geometry.Point([112.56, 26.97])
//var roi = geometry
var roi = ee.FeatureCollection(ee.List([ee.Feature(p1).set('name','p1'),
                                        ee.Feature(p2).set('name','p2')])
                                      )
// Map.addLayer(roi)

//哨兵2号数据
var dataset = ee.ImageCollection("COPERNICUS/S2")
                  .filter(ee.Filter.date('2016-01-01', '2022-10-31'));

1.2.2 去云处理

//去云处理
function maskS2clouds(image) {
  var qa = image.select('QA60')
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;//cirrus clouds:卷云
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0).and(
             qa.bitwiseAnd(cirrusBitMask).eq(0))
  return image.updateMask(mask)
      .select("B.*")
      .copyProperties(image, ["system:time_start"])
}

var filtered = dataset.filter(ee.Filter.bounds(roi))
var filtered = filtered.map(maskS2clouds)
function NDVI(image) {
var ndvi = image.expression(  
    '(NIR-Red) / (NIR+Red)',{  
      'NIR': image.select('B8'),  
      'Red': image.select('B4')
    })
  return image.addBands(ndvi.rename('NDVI')).clip(roi);
}
var withNDVI = filtered.map(NDVI);
print(withNDVI)
var precipitation = withNDVI.select('NDVI');
var ft = ee.FeatureCollection(ee.List([]))
var fill = function(img, ini) {
  var inift = ee.FeatureCollection(ini)
  var ft2 = img.sampleRegions({
  collection:roi,
  properties:ee.List(['name']),
  scale:10
  });
  var date = img.date().format()
  var ft3 = ft2.map(function(f){return f.set("date", date)})
  return inift.merge(ft3)
}

1.2.3 数据获取

// Iterates over the ImageCollection
var newft = ee.FeatureCollection(precipitation.iterate(fill, ft))
Export.table.toDrive({
  collection: newft,
  description: 'sample_get',
  fileFormat: 'CSV'
});

var NDVIchart = ui.Chart.image.series({
                imageCollection: precipitation,
                region:roi,
                reducer: ee.Reducer.mean(),
                scale: 10
              })
              .setOptions({
                title: "NDVI列表", 
                hAxis: {title: "date"},
                vAxis: {title: "ndvi"},
                lineWidth:1,
                pointSize:2
              });
  print("NDVIchart", NDVIchart);

图像展示：

CSV数据：

2. Python多项式拟合

import pandas as pd
import numpy as np
from datetime import datetime
import matplotlib.pyplot as plt
import matplotlib.dates as mdates

data = pd.read_csv('D:/GoogleChrom/points_values.csv')
data = pd.DataFrame(data)
print(data)

date = data.iloc[0:114, 2:3].values.flatten()
ndvi = data.iloc[0:114, 1:2].values.flatten()
date

def d_to_jd(time):
    fmt = '%Y-%m-%d'
    dt = datetime.strptime(time, fmt)
    tt = dt.timetuple()
    return tt.tm_yday
def jd_to_time(time):
    dt = datetime.strptime(time,'%Y%j').date()
    fmt = '%Y/%m/%d'
    return dt.strftime(fmt)

def pol(date,ndvi):
    years = [datetime.strptime(d,'%Y/%m/%d').date() for d in date]
    plt.gca().xaxis.set_major_formatter(mdates.DateFormatter('%m/%d'))
    plt.gca().xaxis.set_major_locator(mdates.DayLocator())
    plt.xticks(years[::30])
    s = list(map(str, years))
    s = list(map(d_to_jd, s))
    _x = np.array(s)
    ndvi = np.array(list(ndvi))
    o = np.polyfit(_x, ndvi, 6)
    g = np.poly1d(o)
#     print(g)
    plt.plot(years, g(_x), c='g', ls='-')
    plt.plot(years, ndvi, c='r', ls='-', marker='o', markersize=2,alpha=0.3)
    plt.legend(labels = ['fitted values','ndvi values'])
    plt.show()

result = pol(date,ndvi)

拟合方程：1.567e-14 x - 1.626e-11 x + 6.464e-09 x - 1.224e-06 x + 0.0001093 x - 0.003585 x + 0.02839

拟合图像结果：