GEE-TED: A tsetse ecological distribution model for Google Earth Engine

GEE-TED: A tsetse ecological distribution model for Google Earth Engine Associated publication forthcoming: Fox, L., Peter, B. G., Frake, A. N., and Messina, J. P. (Forthcoming). A Bayesian maximum entropy model for predicting tsetse ecological distributions. Journal TBD. Description GEE-TED is a Google Earth Engine (GEE; Gorelick et al. 2017) adaptation of a tsetse ecological distribution (TED) model developed by DeVisser et al. (2010), which was designed for use in ESRI's ArcGIS. TED uses time-series climate and land-use/land-cover (LULC) data to predict the probability of tsetse presence across space based on species habitat preferences (in this case Glossina Morsitans). Model parameterization includes (1) day and night temperatures (MODIS Land Surface Temperature; MOD11A2), (2) available moisture/humidity using a vegetation index as a proxry (MODIS NDVI; MOD13Q1), (3) LULC (MODIS Land Cover Type 1; MCD12Q1), (4) year selections, and (5) fly movement rate (meters/16-days). TED has also been used as a basis for the development of an agent-based model by Lin et al. (2015) and in a cost-benefit analysis of tsetse control in Tanzania by Yang et al. (2017). Parameterization in Fox et al. (Forthcoming): Suitable LULC types and climate thresholds used here are specific to Glossina Morsitans in Kenya and are based on the parameterization selections in DeVisser et al. (2010) and DeVisser and Messina (2009). Suitable temperatures range from 17–40°C during the day and 10–40°C at night and available moisture is characterized as NDVI > 0.39. Suitable LULC comprises predominantly woody vegetation; a complete list of suitable categories is available in DeVisser and Messina (2009). In the Fox et al. (Forthcoming) publication, two versions of MCD12Q1 were used to assess suitable LULC types: Versions 051 and 006. The GeoTIFF supplied in this dataset entry (GEE-TED_Kenya_2016-2017.tif) uses the aforementioned parameters to show the probable tsetse distribution across Kenya for the years 2016-2017. A static graphic of this GEE-TED output is shown below and an interactive version can be viewed at: https://cartoscience.users.earthengine.app/view/gee-ted. Figure associated with Fox et al. (Forthcoming) GEE code The code supplied below is generalizable across geographies and species; however, it is highly recommended that parameterization is given considerable attention to produce reliable results. Note that output visualization on-the-fly will take some time and it is recommended that results be exported as an asset within GEE or exported as a GeoTIFF. Note: Since completing the Fox et al. (Forthcoming) manuscript, GEE has removed Version 051 per NASA's deprecation of the product. The current release of GEE-TED now uses only MCD12Q1 Version 006; however, alternative LULC data selections can be used with minimal modification to the code. // Input options var tempMin = 10 // Temperature thresholds in degrees Celsius var tempMax = 40 var ndviMin = 0.39 // NDVI thresholds; proxy for available moisture/humidity var ndviMax = 1 var movement = 500 // Fly movement rate in meters/16-days var startYear = 2008 // The first 2 years will be used for model initialization var endYear = 2019 // Computed probability is based on startYear+2 to endYear var country = 'KE' // Country codes - https://en.wikipedia.org/wiki/List_of_FIPS_country_codes var crs = 'EPSG:32737' // See https://epsg.io/ for appropriate country UTM zone var rescale = 250 // Output spatial resolution var labelSuffix = '02052020' // For file export labeling only //[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17] MODIS/006/MCD12Q1 var lulcOptions006 = [1,1,1,1,1,1,1,1,1, 0, 1, 0, 0, 0, 0, 0, 0] // 1 = suitable 0 = unsuitable // No more input required ------------------------------ // var region = ee.FeatureCollection(\"USDOS/LSIB_SIMPLE/2017\") .filterMetadata('country_co', 'equals', country) // Input parameter modifications var tempMinMod = (tempMin+273.15)/0.02 var tempMaxMod = (tempMax+273.15)/0.02 var ndviMinMod = ndviMin*10000 var ndviMaxMod = ndviMax*10000 var ndviResolution = 250 var movementRate = movement+(ndviResolution/2) // Loading image collections var lst = ee.ImageCollection('MODIS/006/MOD11A2').select('LST_Day_1km', 'LST_Night_1km') .filter(ee.Filter.calendarRange(startYear,endYear,'year')) var ndvi = ee.ImageCollection('MODIS/006/MOD13Q1').select('NDVI') .filter(ee.Filter.calendarRange(startYear,endYear,'year')) var lulc006 = ee.ImageCollection('MODIS/006/MCD12Q1').select('LC_Type1') // Lulc mode and boolean reclassification var lulcMask = lulc006.mode().remap([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17],lulcOptions006) .eq(1).rename('remapped').clip(region) // Merge NDVI and LST image collections var combined = ndvi.combine(lst, true) var combinedList = combined.toList(10000) // Boolean reclassifications (suitable/unsuitable) for day/night temperatures and ndvi var con = function(image) { var ndviExp = image.expression(\"(b('NDVI') > \"+ndviMaxMod+\") ? 0\" + \":(b('NDVI') > \"+ndviMinMod+\") ? 1\" + \":(b('NDVI') < \"+ndviMinMod+\") ? 0\" + \":1\").rename('ndvi_suit') var lstDayExp = image.expression(\"(b('LST_Day_1km') > \"+tempMaxMod+\") ? 0\" + \":(b('LST_Day_1km') > \"+tempMinMod+\") ? 1\" + \":(b('LST_Day_1km') < \"+tempMinMod+\") ? 0\" + \":1\").rename('day_suit') var lstNightExp = image.expression(\"(b('LST_Night_1km') > \"+tempMaxMod+\") ? 0\" + \":(b('LST_Night_1km') > \"+tempMinMod+\") ? 1\" + \":(b('LST_Night_1km') < \"+tempMinMod+\") ? 0\" + \":1\").rename('night_suit') var add = image.addBands(ndviExp).addBands(lstDayExp).addBands(lstNightExp) var multiply = add.select('ndvi_suit').multiply(add.select('day_suit')).multiply(add.select('night_suit')) return add.addBands(multiply.rename('suit')).select('suit') } var conList = combined.map(con).toList(10000) var finish = conList.size().subtract(1).getInfo() // Fly movement rate model var iterateList = conList var kernel = ee.Kernel.square({radius: movementRate, units: 'meters'}) var fill = ee.Image(conList.get(0)).multiply(0) var expansion = fill.add(1) for (var range = 0; range <= finish; range = range + 1) { var img = ee.Image(iterateList.get(range)).select('suit').gt(0).clip(region) var multi = img.multiply(expansion).set('num',range) var mask = multi.eq(1) var masked = multi.updateMask(mask) var expand = masked.focal_max({kernel: kernel}) var expansion = fill.where(expand,1) var updateList = iterateList.add(multi) var iterateList = updateList } // Filter out first two years and compute tsetse probability var filtered = ee.ImageCollection(updateList).sort('num',false).limit(updateList.size().divide(2).subtract(46)) var filteredProb = filtered.sum().divide(filtered.size()).multiply(100).updateMask(lulcMask) var probMask = filteredProb.neq(0) var filteredProbMasked = filteredProb.updateMask(probMask).rename('probability') var startMod = startYear+2 var naming = 'GEE-TED_'+country+'-'+startMod+'-'+endYear+'_'+labelSuffix Map.addLayer(lulcMask, {min:0, max:1}, 'MCD12Q1 LULC mask', false) Map.addLayer(filteredProbMasked, {palette: ['FCFAE1','F6E497','BD8D46','B9121B','4C1B1B'], opacity: 0.8, min: 0, max: 100}, naming, false) Map.setOptions('HYBRID') Map.centerObject(region) // Export image to drive Export.image.toDrive({ image: filteredProbMasked, description: naming, fileNamePrefix: naming, maxPixels: 1e13, scale: rescale, crs: crs, region: region }) References • DeVisser, M.H. and Messina, J.P., 2009. Optimum land cover products for use in a Glossina-morsitans habitat model of Kenya. International Journal of Health Geographics, 8(1), pp.1-20. • DeVisser, M.H., Messina, J.P., Moore, N.J., Lusch, D.P. and Maitima, J., 2010. A dynamic species distribution model of Glossina subgenus Morsitans: The identification of tsetse reservoirs and refugia. Ecosphere, 1(1), pp.1-21. • Gorelick, N., Hancher, M., Dixon, M., Ilyushchenko, S., Thau, D. and Moore, R., 2017. Google Earth Engine: Planetary-scale geospatial analysis for everyone. Remote Sensing of Environment, 202, pp.18-27. • Lin, S., DeVisser, M.H. and Messina, J.P., 2015. An agent-based model to simulate tsetse fly distribution and control techniques: a case study in Nguruman, Kenya. Ecological Modelling, 314, pp.80-89. • Yang, A., Messina, J.P., Grady, S.C. and White, R.A., 2017. Cost–Benefit Analysis of Tsetse Fly Control in Tanzania. Papers in Applied Geography, 3(2), pp.182-195.

### GEE-TED：面向Google Earth Engine的采采蝇生态分布模型
关联待刊论文：Fox, L.、Peter, B. G.、Frake, A. N. 与 Messina, J. P.（待刊）. 《用于预测采采蝇生态分布的贝叶斯最大熵模型》，期刊待定。

## 数据集概述
GEE-TED是Google Earth Engine（GEE; Gorelick et al., 2017）对DeVisser等人（2010）开发的采采蝇生态分布（Tsetse Ecological Distribution, TED）模型的适配版本，该原始模型最初设计用于ESRI的ArcGIS平台。TED模型通过时序气候与土地利用/土地覆盖（Land Use/Land Cover, LULC）数据，基于物种栖息地偏好（本研究针对舌蝇（Glossina Morsitans））预测空间范围内的采采蝇存在概率。

模型参数化内容包括：（1）昼夜地表气温（采用MODIS地表温度产品MOD11A2）；（2）以植被指数为代理变量的可用湿度/水汽条件（采用MODIS归一化植被指数NDVI产品MOD13Q1）；（3）土地利用/土地覆盖数据（采用MODIS土地覆盖类型1产品MCD12Q1）；（4）研究年份范围选择；（5）蝇类移动速率（单位：米/16天）。

TED模型曾被Lin等人（2015）用作基于智能体的模型开发的基础框架，也被Yang等人（2017）用于坦桑尼亚采采蝇防控的成本效益分析。

### Fox等人（待刊）的参数设置
本研究采用的适宜LULC类型与气候阈值针对肯尼亚境内的舌蝇（Glossina Morsitans），基于DeVisser等人（2010）与DeVisser和Messina（2009）的参数选择方案。适宜气温范围为日间17~40℃，夜间10~40℃；可用湿度以NDVI>0.39作为表征指标。适宜的LULC类型以木本植被为主，完整的适宜类别列表可参见DeVisser与Messina（2009）。在Fox等人（待刊）的研究中，共使用了两个版本的MCD12Q1数据以评估适宜LULC类型：版本051与版本006。

本数据集条目内提供的GeoTIFF文件（`GEE-TED_Kenya_2016-2017.tif`）采用上述参数，展示了2016-2017年肯尼亚境内的潜在采采蝇分布。该GEE-TED输出结果的静态可视化图见下文，交互式版本可访问：https://cartoscience.users.earthengine.app/view/gee-ted。该图为Fox等人（待刊）研究的配套附图。

## GEE代码说明
下文提供的代码可在不同地理区域与物种间通用，但强烈建议仔细调整参数以获得可靠的分析结果。请注意，实时可视化输出需要一定计算时间，建议将结果导出为GEE内置资产或GeoTIFF格式。

注意：在完成Fox等人（待刊）的手稿后，GEE已根据NASA对该产品的停用声明移除了MCD12Q1版本051。当前发布的GEE-TED仅使用MCD12Q1版本006；不过，只需对代码进行少量修改即可接入其他LULC数据集。

javascript
// 输入参数设置
var tempMin = 10 // 以摄氏度为单位的最低适宜气温
var tempMax = 40 // 以摄氏度为单位的最高适宜气温
var ndviMin = 0.39 // NDVI阈值，作为可用湿度/水汽的代理变量
var ndviMax = 1
var movement = 500 // 蝇类移动速率，单位：米/16天
var startYear = 2008 // 初始前2年用于模型初始化
var endYear = 2019 // 计算得到的概率基于startYear+2至endYear的数据
var country = 'KE' // 国家代码，详见：https://en.wikipedia.org/wiki/List_of_FIPS_country_codes
var crs = 'EPSG:32737' // 坐标系，详见https://epsg.io/获取对应国家的UTM分区
var rescale = 250 // 输出空间分辨率
var labelSuffix = '02052020' // 仅用于文件导出的命名后缀
// MODIS/006/MCD12Q1的类别索引：[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]
var lulcOptions006 = [1,1,1,1,1,1,1,1,1, 0, 1, 0, 0, 0, 0, 0, 0] // 1=适宜，0=不适宜
// 无需额外输入 ------------------------------

// var region = ee.FeatureCollection("USDOS/LSIB_SIMPLE/2017")
// .filterMetadata('country_co', 'equals', country)

// 输入参数修正
var tempMinMod = (tempMin+273.15)/0.02
var tempMaxMod = (tempMax+273.15)/0.02
var ndviMinMod = ndviMin*10000
var ndviMaxMod = ndviMax*10000
var ndviResolution = 250
var movementRate = movement+(ndviResolution/2)

// 加载影像集合
var lst = ee.ImageCollection('MODIS/006/MOD11A2').select('LST_Day_1km', 'LST_Night_1km')
.filter(ee.Filter.calendarRange(startYear,endYear,'year'))
var ndvi = ee.ImageCollection('MODIS/006/MOD13Q1').select('NDVI')
.filter(ee.Filter.calendarRange(startYear,endYear,'year'))
var lulc006 = ee.ImageCollection('MODIS/006/MCD12Q1').select('LC_Type1')

// 土地覆盖模式与布尔重分类
var lulcMask = lulc006.mode().remap([1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17],lulcOptions006)
.eq(1).rename('remapped').clip(region)

// 合并NDVI与LST影像集合
var combined = ndvi.combine(lst, true)
var combinedList = combined.toList(10000)

// 对昼夜气温、NDVI进行布尔重分类（适宜/不适宜）
var con = function(image) {
var ndviExp = image.expression("(b('NDVI') > "+ndviMaxMod+") ? 0" + ":(b('NDVI') > "+ndviMinMod+") ? 1" + ":(b('NDVI') < "+ndviMinMod+") ? 0" + ":1").rename('ndvi_suit')
var lstDayExp = image.expression("(b('LST_Day_1km') > "+tempMaxMod+") ? 0" + ":(b('LST_Day_1km') > "+tempMinMod+") ? 1" + ":(b('LST_Day_1km') < "+tempMinMod+") ? 0" + ":1").rename('day_suit')
var lstNightExp = image.expression("(b('LST_Night_1km') > "+tempMaxMod+") ? 0" + ":(b('LST_Night_1km') > "+tempMinMod+") ? 1" + ":(b('LST_Night_1km') < "+tempMinMod+") ? 0" + ":1").rename('night_suit')
var add = image.addBands(ndviExp).addBands(lstDayExp).addBands(lstNightExp)
var multiply = add.select('ndvi_suit').multiply(add.select('day_suit')).multiply(add.select('night_suit'))
return add.addBands(multiply.rename('suit')).select('suit')
}
var conList = combined.map(con).toList(10000)
var finish = conList.size().subtract(1).getInfo()

// 蝇类移动速率模型
var iterateList = conList
var kernel = ee.Kernel.square({radius: movementRate, units: 'meters'})
var fill = ee.Image(conList.get(0)).multiply(0)
var expansion = fill.add(1)
for (var range = 0; range <= finish; range = range + 1) {
var img = ee.Image(iterateList.get(range)).select('suit').gt(0).clip(region)
var multi = img.multiply(expansion).set('num',range)
var mask = multi.eq(1)
var masked = multi.updateMask(mask)
var expand = masked.focal_max({kernel: kernel})
var expansion = fill.where(expand,1)
var updateList = iterateList.add(multi)
var iterateList = updateList
}

// 过滤前两年数据并计算采采蝇存在概率
var filtered = ee.ImageCollection(updateList).sort('num',false).limit(updateList.size().divide(2).subtract(46))
var filteredProb = filtered.sum().divide(filtered.size()).multiply(100).updateMask(lulcMask)
var probMask = filteredProb.neq(0)
var filteredProbMasked = filteredProb.updateMask(probMask).rename('probability')
var startMod = startYear+2
var naming = 'GEE-TED_'+country+'-'+startMod+'-'+endYear+'_'+labelSuffix

Map.addLayer(lulcMask, {min:0, max:1}, 'MCD12Q1 LULC掩膜', false)
Map.addLayer(filteredProbMasked, {palette: ['FCFAE1','F6E497','BD8D46','B9121B','4C1B1B'], opacity: 0.8, min: 0, max: 100}, naming, false)
Map.setOptions('HYBRID')
Map.centerObject(region)

// 将影像导出至Google云端硬盘
Export.image.toDrive({
image: filteredProbMasked,
description: naming,
fileNamePrefix: naming,
maxPixels: 1e13,
scale: rescale,
crs: crs,
region: region
})

## 参考文献
1. DeVisser, M.H. 与 Messina, J.P., 2009. 适用于肯尼亚舌蝇（Glossina Morsitans）栖息地模型的最优土地覆盖产品. 《国际卫生地理学杂志》, 8(1), 第1-20页.
2. DeVisser, M.H.、Messina, J.P.、Moore, N.J.、Lusch, D.P. 与 Maitima, J., 2010. 舌蝇亚属（Morsitans）的动态物种分布模型：采采蝇储源地与避难所的识别. 《生态圈》, 1(1), 第1-21页.
3. Gorelick, N.、Hancher, M.、Dixon, M.、Ilyushchenko, S.、Thau, D. 与 Moore, R., 2017. Google Earth Engine：面向大众的行星级地理空间分析平台. 《环境遥感》, 202, 第18-27页.
4. Lin, S.、DeVisser, M.H. 与 Messina, J.P., 2015. 模拟采采蝇分布与防控技术的基于智能体的模型：以肯尼亚恩古鲁曼为例. 《生态建模》, 314, 第80-89页.
5. Yang, A.、Messina, J.P.、Grady, S.C. 与 White, R.A., 2017. 坦桑尼亚采采蝇防控的成本效益分析. 《应用地理学论文集》, 3(2), 第182-195页.