看小伙是如何用python尽览我国植被覆盖率的沧海变化

前言

昨天学弟给了我一个数据包“我国NDVI的栅格数据”（即归一化植被指数），需要把陕西省的数据批量裁剪出来。

于是，我拿出python。“刷刷刷~”，搞定了：

这就完事了？ 作为深度践行勤俭节约的社会主义好青年的kimol君，当然不会浪费资源。于是，我决定用这个数据来分析一下我国植被覆盖率变化的沧海桑田（以陕西省为例）。

先上图，让大家伙一睹为快：

可以很直观地看出，从1998-2018年，陕西省（尤其是陕北地区）的植被覆盖率正在不断增加。我想，这离不开国家治沙工程的大力推行，更离不开许许多多护林治沙先辈们的辛勤付出~（手动敬礼中…）

一、栅格数据转为灰度图

本次用到的数据是NDVI栅格数据，它是以“.tif”为后缀的。简单来理解的话，我们可以把它看作一张图，每个像素点的值代表了该位置的归一化植被指数。为了更方便展示，我们将其转为灰度图：

首先读取tif数据：

def transfer(img):
    '''
    将tif转为灰度图
    '''
    for x in range(img.shape[0]):
        for y in range(img.shape[1]):
            if img[x][y] < 0:
                img[x][y] = 0
            else:
                img[x][y] *= 255
    return img
    
# 创建文件夹用于存储灰度图
grayPath = './gray'
if not os.path.exists(grayPath):
    os.mkdir(grayPath)

# 遍历tif文件
tifPath ='./data'
for file in os.listdir(tifPath):
    if file.endswith('.tif'):        
        img = cv2.imread('%s/%s'%(tifPath, file), 2)
        transfer(img) # 转为灰度格式
        cv2.putText(img,'Year %s'%file[9:13], (20,40), cv2.FONT_HERSHEY_SIMPLEX, 1.0,(255,255,255), 1, cv2.LINE_AA) # 加入年份说明
        cv2.imwrite('%s/%s.png'%(grayPath, file.split('.')[0]), img)
        print('"%s"处理完毕!'%file)

其中cv2.putText是为了在图片上加入年份，方便区分，最后得到灰度图如下：

然而，灰度图似乎显得不太友好，为了更直观地表现出植被覆盖率，我决定将其转为“伪彩色”。

二、灰度图转为伪彩色图

我们根据像素点的灰度值来定义三个颜色区间（红、黄、绿），并进行转换：

def trans_color(img):
    '''
    将灰度图转为伪彩色
    '''
    for x in range(img.shape[0]):
        for y in range(img.shape[1]):
            if img[x][y][0] == 0 and img[x][y][1] == 0 and img[x][y][2] == 0: # 边缘区域跳过
                continue
            gray = img[x][y][0]
            if gray <= 110: # 红
                img[x][y][0] = 0
                img[x][y][1] = gray*1.5
                img[x][y][2] = 255
            elif gray <= 200: # 黄
                gray -= 85
                img[x][y][0] = 0
                img[x][y][1] = 255
                img[x][y][2] = 255-(127*gray/85)
            else: # 绿
                gray -= 170
                img[x][y][0] = 0
                img[x][y][1] = 255
                img[x][y][2] = 127-(127*gray/85)
    return img

需要注意的是，cv2.imread默认读取的模式并非RGB而是BRG。

遍历灰度图片，并转换：

# 创建文件夹用于存储彩色图
rgbPath = './rgb'
if not os.path.exists(rgbPath):
    os.mkdir(rgbPath)
    
# 将灰度图转为彩色图
for file in os.listdir(grayPath):
    img = cv2.imread('%s/%s'%(grayPath, file))
    trans_color(img) # 转换
    cv2.imwrite('%s/%s'%(rgbPath, file), img)
    print('"%s"转换完毕!'%file)

于是，我们得到了彩色的图片：

三、合成为gif

利用imageio库，将众多图片合成为一张gif图：

def create_gif(image_list, gif_name, duration=0.35):
    '''
    生成gif
    '''
    frames = []
    for image_name in image_list:
        frames.append(imageio.imread(image_name))
    imageio.mimsave(gif_name, frames, 'GIF', duration=duration)
    return

# 合成gif图
images_list = []
for file in os.listdir(rgbPath):
    images_list.append('%s/%s'%(rgbPath, file))
create_gif(images_list, 'NDVI.gif')

于是，我们便有了前言中的那张图，收工~

写在最后

先辈们用坚毅的品质告诉了我们什么叫敢为人先；
先辈们用辛勤的汗水告诉了我们什么叫持之以恒；
先辈们用喜人的成果告诉了我们什么叫锐意进取。

我辈又有何种理由不奋进呢？

我是kimol君，咋们下次再会~

创作不易，大侠请留步… 动起可爱的双手，来个赞再走呗 (๑◕ܫ￩๑)