作物模型的区域模拟实现过程初探

作物生长模型能够根据作物的生长环境模拟作物的生长发育，为农业生产及科学研究提供有价值的参考。但模型多数是基于单点或者几个站点的模拟结果进行差值来表征区域的模拟结果。这里来探讨一下模型的网格化模拟实现。
一般的作物生长模型都需要土壤、气象、管理、初始条件等几个必要的因素。在区域模拟时还需要种植区域的栅格图。

基本流程

假设我们现在有一副小麦区域分布图 wheat.tif
以下以Python代码为例提供一些过程的核心代码

种植区域获取

首先需要获取种植区域的坐标
将栅格的种植区域读取为array

"""
author: Shuai-jie Shen 沈帅杰
CSDN: https://blog.csdn.net/weixin_45452300
公众号: AgBioIT
"""
from osgeo import gdal
import numpy
import  pandas  as pd


gdal.AllRegister()
filePath = r'wheat.tif'
dataset = gdal.Open(filePath)
adfGeoTransform = dataset.GetGeoTransform()
nXSize = dataset.RasterXSize  # 列数
nYSize = dataset.RasterYSize  # 行数
im_data = dataset.ReadAsArray(0,0,nXSize,nYSize)
index = []  # 纬度
columns = []  # 经度
for j in range(nYSize):
    lat = adfGeoTransform[3] + j * adfGeoTransform[5]
    index.append(lat)
for i in range(nXSize):
    lon = adfGeoTransform[0] + i * adfGeoTransform[1]
    columns.append(lon)
data = pd.DataFrame(im_data, index=index, columns=columns)

此时我们得到了栅格化的种植区域图dataframe，假设种植区域值为1，非种植区域为0
然后对结果进行遍历，找到为1的点根据根据坐标进行模型模拟，为0 的直接跳过（ps:这个方法应该是比较笨的，但目前我认知只能实现这个了）

# 创建一个空列表储存模拟结果
results = pd.DataFrame(np.zeros(nYSize, nXSize), index=index, columns=columns)
for i in data.index:
	for j in data.columns:
		if data.loc[i, j] == 0:
			continue
		runmodel  # 运行自己的模型
		results .loc[i, j] = runmodel.result # 写入运行结果

根据栅格坐标run model

气象数据选择

气象数据可以直接从NASA下载，NASA提供了1984年至今全球0.5度分辨率的气象数据
这个可以提前下载好覆盖整个区域的数据，可以参考从NASA获取全球气象数据
这里提供一个函数，能够将坐标值近似到0.5度的函数

def st_loc(num):
    """
    选取气象数据的时候标准化输入为0.5的分辨率
    """
    import math
    if num-math.floor(num) <=0.25 or num-math.floor(num)>=0.75:
        result = round(num)
    elif 0.25 < num-math.floor(num) < 0.75:
        result = math.floor(num) + 0.5      
    return result

土壤数据的选择

土壤数据可以从数据库中获取（这里选取青藏高原数据科学中心的中国数据集），此数据集提供了1km和250m分辨率的全国土壤特征，土层深度有0-20, 0-5, 5-15, 15-30, 30-60, 60-100, 100-200cm，但也需要一定的转化，土壤水分参数估计可以参考土壤水分特征参数估计
土壤数据为tif文件格式，读取的时候与种植区域获取的一样，直接读取为array，index和columns分别是经纬度。
但这个经纬度和种植区域的经纬度可能不完全一致，因此可以利用以下函数，检索最近的位置。利用此函数就可以根据种植图的坐标定位到土壤数据集的相应坐标位置

from bisect import bisect_left
def takeClosest(myList, myNumber):
    if myNumber >= max(myList):
        return max(myList)
    elif myNumber <= min(myList):
        return min(myList)
    if myList[0] < myList[1]:
        pos = bisect_left(myList, myNumber)   # 找到 mylist 里面第一个不比 mynumber 小（即 >= ）的数的索引下标
        # 返回的插入点 pos 可以将数组myList分成两部分。左侧是 all(val < x for val in myList[lo:i]) ，右侧是 all(val >= x for val in myList[i:hi])
        before = myList[pos - 1]
        after = myList[pos]
        if after - myNumber < myNumber - before:
           return after
        else:
           return before
    else:
        myList = myList[::-1]
        pos = takeClosest(myList, myNumber)
        return pos

效果

list=[35.11225,35.2365556,35.3656545,35.4556685]
answer = takeClosest(cityarea_list, 35.40)
print(answer)

输出：35.3656545

一些管理措施如灌溉，播种等均可找到响应的数据集，例如：灌溉量数据集当然，模型也可以自己设置管理情景。