Preliminary study on Regional Simulation of crop model

Crop growth model can simulate the growth and development of crops according to the growth environment of crops , It provides valuable reference for agricultural production and scientific research . However, most of the models are based on the difference between the simulation results of a single point or several stations to represent the regional simulation results . Here, let's discuss the grid simulation implementation of the model .
General crop growth models require soil 、 meteorological 、 management 、 Initial conditions And so on . In the area simulation, the grid map of planting area is also required .

The basic flow

Suppose we now have a regional distribution map of wheat wheat.tif
The following Python Code as an example to provide some core code of the process

Planting area acquisition

First, you need to get the coordinates of the planting area
Read the planting area of the grid as array

""" author: Shuai-jie Shen  Shenshuaijie  CSDN: https://blog.csdn.net/weixin_45452300  official account : 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  #  Number of columns 
nYSize = dataset.RasterYSize  #  Row number 
im_data = dataset.ReadAsArray(0,0,nXSize,nYSize)
index = []  #  latitude 
columns = []  #  longitude 
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)

At this point, we get the grid planting area map dataframe, Suppose the planting area value is 1, The non planting area is 0
Then iterate over the result , Find as 1 The points of the model are simulated according to the coordinates , by 0 We'll skip it （ps: This method should be stupid , But now I can only realize this ）

#  Create an empty list to store simulation results 
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  #  Run your own model 
		results .loc[i, j] = runmodel.result #  Write run results 

According to the grid coordinates run model

Meteorological data selection

Meteorological data can be obtained directly from NASA download ,NASA Provides 1984 Year to date 0.5 Degree resolution meteorological data
You can download the data covering the whole area in advance , You can refer to from NASA Access to global meteorological data
Here is a function , Be able to approximate the coordinate value to 0.5 Function of degree

def st_loc(num):
    """  When selecting meteorological data, the standardized input is 0.5 The resolution of the  """
    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

Selection of soil data

Soil data can be obtained from the database （ Here, the Chinese data set of the Qinghai Tibet Plateau Data Science Center is selected ）, This dataset provides 1km and 250m Resolution of national soil characteristics , The depth of soil layer is 0-20, 0-5, 5-15, 15-30, 30-60, 60-100, 100-200cm, But it also needs some transformation , The estimation of soil moisture parameters can refer to Estimation of soil moisture characteristic parameters
The soil data is tif File format , The reading time is the same as that of the planting area , Read directly as array,index and columns They are longitude and latitude .
However, this longitude and latitude may not be completely consistent with the longitude and latitude of the planting area , So you can use the following functions , Retrieve the most recent location . This function can be used to locate the corresponding coordinate position of the soil data set according to the coordinates of the planting map

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)   #  find  mylist  The first one is no better than  mynumber  Small （ namely  >= ） Index subscript of the number of 
        #  The insertion point returned  pos  You can put arrays myList In two parts . On the left is  all(val < x for val in myList[lo:i]) , The right side is  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

effect

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

Output ：35.3656545

Some management measures such as irrigation , The response data set can be found by seeding , for example ： Irrigation volume data set Of course , The model can also set its own management scenarios .