Modis-ET-NPP-GPP

Name: modis数据处理

Description: modis数据读取、去除异常值、投影转换、裁剪和栅格求和（蒸散发数据、NPP、GPP）

Parameters 参数1: modis影像 #input,infile,.hdf 参数2: 输出影像 #output,outfile,.tif

read_img函数：读栅格文件，返回影像的im_height, im_width, im_geotrans, im_proj, im_data；

write_img函数：保存栅格文件；

readHdfWithGeo函数：读取hdf中的子数据集，并去除异常值，保存为tif；

cut_img函数：按照矢量进行裁剪；

sum_img函数：8天一期的栅格数据求和，得到一年的累积值。

#!/usr/bin/env python3
# -*- coding: utf-8 -*-"""
Name:modis数据处理
Description:modis数据读取、去除异常值、投影转换、裁剪和栅格求和（蒸散发数据、NPP、GPP）
Parameters参数1: modis影像 #input,infile,*.hdf参数2: 输出影像 #output,outfile,*.tif
"""import os
import os.path
import sys
from osgeo import gdal, ogr
from osgeo import gdal_array
import time
import numpy as np# ignore warning
np.seterr(divide='ignore', invalid='ignore')
gdal.SetConfigOption('CPL_LOG', 'NUL')
gdal.UseExceptions()
ogr.UseExceptions()start = time.time()
driver_img = gdal.GetDriverByName('GTiff')
driver_shp = ogr.GetDriverByName('ESRI Shapefile')# 读栅格文件
def read_img(filename):dataset = gdal.Open(filename)  # 打开文件print(dataset)im_width = dataset.RasterXSize  # 栅格矩阵的列数print(im_width)im_height = dataset.RasterYSize  # 栅格矩阵的行数im_bands = dataset.RasterCount  # 波段数im_geotrans = dataset.GetGeoTransform()  # 仿射矩阵，左上角像素的大地坐标和像素分辨率im_proj = dataset.GetProjection()  # 地图投影信息，字符串表示im_data = dataset.ReadAsArray(0, 0, im_width, im_height)  # 将数据写成数组，对应栅格矩阵del datasetreturn im_height, im_width, im_geotrans, im_proj, im_data# 写栅格文件
def write_img(filename, im_proj, im_geotrans, im_data):# 判断栅格数据的数据类型if 'int8' in im_data.dtype.name:datatype = gdal.GDT_Byteelif 'int16' in im_data.dtype.name:datatype = gdal.GDT_UInt16else:datatype = gdal.GDT_Float32# 判断数组维数if len(im_data.shape) == 3:im_bands, im_height, im_width = im_data.shapeelse:im_bands, (im_height, im_width) = 1, im_data.shape# 创建文件driver = gdal.GetDriverByName('GTiff')  # 数据类型必须有dataset = driver.Create(filename, im_width, im_height, im_bands, datatype)dataset.SetGeoTransform(im_geotrans)  # 写入仿射变换参数dataset.SetProjection(im_proj)  # 写入投影if im_bands == 1:dataset.GetRasterBand(1).WriteArray(im_data)else:for i in range(im_bands):dataset.GetRasterBand(i + 1).WriteArray(im_data[i])del dataset# 读取hdf中的波段，保存为tif
def readHdfWithGeo(hdfFloder, saveFloder, refer_data_file):hdfNameList = os.listdir(hdfFloder)print('文件个数：', len(hdfNameList))refer_data = gdal.Open(refer_data_file)proj = refer_data.GetProjectionRef()im_geotrans = refer_data.GetGeoTransform()# 遍历文件名列表中所有文件for i in range(len(hdfNameList)):dirname, basename = os.path.split(hdfNameList[i])  # 获取文件名后缀filename, txt = os.path.splitext(basename)hdfPath = hdfFloder + os.sep + hdfNameList[i]  # os.sep: 跨平台的文件路径分隔符datasets = gdal.Open(hdfPath)sus = datasets.GetSubDatasets()et_rows, et_columns, et_geotrans, et_proj, et_array = read_img(sus[1][0])print('产品子数据集个数为：', len(sus))et_array = et_array.astype(np.float32)et_array[np.where((et_array < -32761) | (et_array > 32761))] = 0  # <-32767以及>32700设为空值，不参与计算Scale_Factor = 0.0001et_array = et_array * Scale_Factorwrite_img(saveFloder + os.sep + filename + '_pro.tif',proj, im_geotrans, et_array)# 栅格裁剪
def cut_img(tifFloder, saveFloder, shpFile):tifNameList = os.listdir(tifFloder)for i in range(len(tifNameList)):filename, txt = os.path.splitext(tifNameList[i])if txt == '.tif':tifPath = tifFloder + os.sep + tifNameList[i]outName = saveFloder + os.sep + filename + '_cut' + txtcut_img = gdal.Warp(outName, tifPath,cutlineDSName=shpFile,cropToCutline=True)del cut_imgprint('{filename} deal end '.format(filename=outName))# 栅格批量求和
def sum_img(tifFloder, saveFloder, refer_data_file):rows, cols, geotrans, proj, data_array = read_img(refer_data_file)filenames = os.listdir(tifFloder)nd_arr = np.zeros((len(filenames), rows, cols))i = 0for filename in filenames:if os.path.splitext(filename)[1] == '.tif':dataset = gdal.Open(tifFloder + os.sep + filename)array = dataset.ReadAsArray()nd_arr[i] = arraydel arrayi += 1sum_arr = np.sum(nd_arr, axis=0)write_img(saveFloder + os.sep + filename + '_sum.tif',proj, geotrans, sum_arr)if __name__ == "__main__":start = time.time()# 参考影像和矢量refer_data_file = r'D:\Data\data_modis\data_results\MOD17A2\MOD17A2_clip\MOD17A2H.A2021361.h27v05.006.2022005050243_pro_cut.tif'shpFile = r'D:\Data\data_modis\data_reference\hebi_pro.shp'# 文件路径hdfFloder = r'D:\Data\data_modis\data_ori\MOD173AHGFv061'proFloder = r'D:\Data\DEM\dem_500m'clipFloder = r'D:\Data\data_modis\data_results\MOD17A2\MOD17A2_clip'saveFloder = r'D:\Data\data_modis\data_results\MOD17A2\sum_img'# 执行函数# readHdfWithGeo(hdfFloder, proFloder, refer_data_file)# cut_img(proFloder, clipFloder, shpFile)sum_img(clipFloder, saveFloder, refer_data_file)end = time.time()print('deal spend: {s} s'.format(s=end - start))

示例数据下载地址：

cll011/Modis-ET-NPP-GPP (github.com)

Modis-ET-NPP-GPP相关推荐

GIS实战应用案例100篇（二十一）-全国分省、市、县净初级生产力NPP数据制作实战（附代码）
前言净初级生产力(Net primary productivity, NPP)是研究陆地生态系统中物质和能量转换的重要指标,NPP的空间分布与区域气候.植被生长以及人类活动等因素息息相关,其变化能反 ...
asp 取数据计算_地学数据 | 地理空间数据获取方式汇总
1.测绘地理信息局会 (http://www.webmap.cn/main.do?method=index) 该网站提供:30米全球地表覆盖数据,GlobeLand30能够提供包括:地理位置.分布范围 ...
转|地理数据下载的网站汇总全国地研联测绘学报 2018-12-05
国内 http://www.dsac.cn/ 地理国情监测云平台 http://henu.geodata.cn/index.jsp 黄河流域数据中心 http://www.loess.csdb.cn/ ...
地理数据下载网站汇总
转自:微信公众号"测绘学报" 国内 http://www.dsac.cn/ 地理国情监测云平台 http://henu.geodata.cn/index.jsp 黄河流域数据中心 ...
Coding and Paper Letter（五十）
资源整理.春节后第一次更文. 文章目录 1 Coding: 2 Paper: 1 Coding: 1.Cartographer是一个跨多个平台和传感器配置提供2D和3D实时同步定位和映射(SLAM)的 ...
matlab 线性回归参数显著性,基于Matlab的栅格数据一元线性回归及显著性检验（slope趋势分析）...
%by yinlichang3064@163.com 在进行长时间序列的栅格数据分析时,如NDVI,fvc,LAI,NPP,GPP,需要知道每个格点的长期趋势. 如果再arcgis中进行一元回归计算, ...
matlab计算栅格数据逐像元趋势,基于Matlab的栅格数据一元线性回归及显著性检验（slope趋势分析）...
%by yinlichang3064@163.com 在进行长时间序列的栅格数据分析时,如NDVI,fvc,LAI,NPP,GPP,需要知道每个格点的长期趋势. 如果再arcgis中进行一元回归计算, ...
Google Earth Engine ——用GPP和NPP计算生物量（朝鲜为例）
数据: MOD17A3H V6产品提供500米像素分辨率的年度净初级生产力(NPP)信息.年净初级生产力是由给定年份的45个8天净光合作用(PSN)产品(MOD17A2H)之和得出的.PSN值是GPP ...
【GEE】基于MODIS产品的NPP NDVI EVI数据提取
其中植被指数的获取,采用MOD13Q1.MOD13A1产品,NPP获取采用MOD17A3HGF产品. MOD13Q1:全称为MODIS/Terra Vegetation Indices 16-Day ...
MODIS影像获取NPP均值
MODIS影像获取NPP均值 var roi = ee.FeatureCollection("users/tnorby415/GD"); var data = ee.ImageCo ...

Modis-ET-NPP-GPP

Modis-ET-NPP-GPP相关推荐

最新文章

热门文章