首页 > 美文鉴赏

python3计算各类地物NDVI均值及将NDVI异常值栅格像元写为点矢量图形

更新时间:2023-07-13 00:58:20 阅读：评论：0

python3计算各类地物NDVI均值及将NDVI异常值栅格像元写为点⽮量图形

请利⽤下列数据统计各类⼟地利⽤类型的NDVI均值。NDVI值通常位于-1⾄1之间，请将该地区NDVI异常值发⽣的像元位置和具体数值，⽣成⼀个点shape⽂件。

数据：

(1) Landsat8地⾯反射率⽂件

LC08_L1TP_121041_121042_20190312_T1_sr_scale_mosaic_clip.tif

Landsat8的NDVI公式：NDVI = (Band5-Band4)/(Band5+Band4)

(2) 梅川江流域的⼟地利⽤图及其说明⽂件

landu30_clip.tif

landu_lookup.csv

考察要点：

栅格⽂件的读写、⽮量shape⽂件的读写（⽐须掌握）

Numpy的熟练应⽤，巧⽤meshgrid函数（⽐须掌握）

Matplotlib绘图（本次作业可选⽤）

提交要求：

（1）基本的python代码⽂件，（2）在ArcMap中，以⼟地利⽤图为背景、以SHAPE点为前景（NDVI正的异常值点为红⾊，负的异常值点为蓝⾊）做⼀幅图，拷屏（或导出map）⽣成⼀个jpg⽂件。如果⼤家积极主动，可以⽤matplotlib把⼟地利⽤图和NDVI异常值shape⽂件叠置在⼀起，⽣成⼀个jpg⽂件。

（3）将上述两个⽂件压缩到⼀个rar⽂件中，命名为"学号_姓名.rar"提交到课程⽹站，如：2019001002_zhangwentian.rar。

请⼤家独⽴完成，作为结课成绩的⼀部分。

说明

landu_lookup. csv 存储了“⼟地利⽤编码”和“SWAT内植被类型”的查找表。其中，

AGRL 农⽥

SWRN 草地

RNGB 灌丛

WATR ⽔体

URHD 建筑⽤地

FRSD 落叶林

FRSE 常绿林

FRST 混交林

landu_code swat-vegname

10 AGRL

11 AGRL

12 AGRL

30 SWRN

40 RNGB致家人的一封信

60 WATR

61 WATR

63 WATR

80 URHD

90 URHD

211 FRSD

221 FRSE

230 FRST

敬天爱人

代码（python3）：

# -*- coding: utf-8 -*-

"""

Created on Sat Oct 19 23:27:54 2019

@author: Depei Bai

"""

from osgeo import gdal,ogr,osr

import numpy as np

import os

import osgeo

#选中代码，tab统⼀右进四格，shift+tab统⼀左退两格

#--------------------------------------read data------------------------------------

#读取⼟地利⽤数据

datat_1 = gdal.Open("landu30_clip.tif")

d1_samples = datat_1.RasterXSize

d1_lines = datat_1.RasterYSize

d1_bands = datat_1.RasterCount

geotans_1 = datat_1.GetGeoTransform()

proj_1 = datat_1.GetProjection()

data_1 = datat_1.ReadAsArray(0,0,d1_samples,d1_lines)

del datat_1

#读取OLI数据，计算NDVI图像

datat_2 = gdal.Open("LC08_L1TP_121041_121042_20190312_T1_sr_scale_mosaic_clip.tif")

d2_samples = datat_2.RasterXSize

d2_lines = datat_2.RasterYSize

d2_bands = datat_2.RasterCount

geotans_2 = datat_2.GetGeoTransform()#左上⾓像元的⼤地坐标和图像分辨率

# 如果图像不含地理坐标信息，默认返回值是：(0,1,0,0,0,1)

#In a north up image,

#padfGeoTransform即为datat.GetGeoTransform()返回的值

#padfGeoTransform[0],padfGeoTransform[3]指左上⾓点坐标(x，y)；

#padfGeoTransform[1]图像在x⽅向的空间分辨率，padfGeoTransform[5]是图像在y⽅向上的空间分辨率；

#padfGeoTransform[2]和padfGeoTransform[4]指正北⽅向旋转⾓度

#如果图像正上指北的，则这两个参数为0

#Fetches the coefficients for transforming between 列/⾏ (P,L) raster space, and projection coordinates (Xp,Yp) space #Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];

#Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];

#————————————————

绩效诊断

proj_2 = datat_2.GetProjection()#指椭球体和投影⽅式

data_2 = datat_2.ReadAsArray(0,0,d1_samples,d1_lines)#取了所有波段所有⾏列的数据

del datat_2

#----------------------------------Calculate NDVI-------------------------------

NDVI =(data_2[4,:,:]*1.0- data_2[3,:,:]*1.0)/(data_2[4,:,:]+ data_2[3,:,:])

#----------------------------------calculate the meam_NDVI_value---------------

#每种地物类型统计其NDVI值的时候均需要取⾮异常值进⾏统计，剔除该类型中⼤于1和⼩于0的值

#np.where()⽆法在⼀个括号中单纯添加多个条件，⽽是要⽤多个括号，& 或 |

#APPEND⽆返回值，因此不能⽤赋值符号

#list to numpy.array强制转换时不能在原值的基础上改，只能重新给到新的变量

#⽤作索引的数组必须是整型或者布尔型，np.array型是不能作为下标的，元组也不可以

AGRL = np.where(((data_1 ==10)|(data_1 ==11)|(data_1 ==12))&(NDVI >-1)&(NDVI <1))

mean_agrl = np.mean(NDVI[AGRL])

SWRN = np.where((data_1 ==30)&(NDVI >-1)&(NDVI <1))

mean_swrn = np.mean(NDVI[SWRN])

RNGB = np.where((data_1 ==40)&(NDVI >-1)&(NDVI <1))

mean_rngb = np.mean(NDVI[RNGB])

WATR = np.where(((data_1 ==60)|(data_1 ==61)|(data_1 ==63))&(NDVI >-1)&(NDVI <1))

mean_watr = np.mean(NDVI[WATR])

URHD = np.where(((data_1 ==80)|(data_1 ==90))&(NDVI >-1)&(NDVI <1))

mean_urhd = np.mean(NDVI[URHD])

FRSD = np.where((data_1 ==211)&(NDVI >-1)&(NDVI <1))

mean_frsd = np.mean(NDVI[FRSD])

FRSE = np.where((data_1 ==221)&(NDVI >-1)&(NDVI <1))

mean_fr = np.mean(NDVI[FRSE])

FRST = np.where((data_1 ==230)&(NDVI >-1)&(NDVI <1))

mean_frst = np.mean(NDVI[FRST])

positive_outlier_index = np.where(NDVI >1)

negative_outlier_index = np.where(NDVI <-1)

#---------------------------将异常值写⼊shape⽂件中------------------------------

osgeo.gdal.SetConfigOption('GDAL_FILENAME_IS_UTF8','NO')#解决中⽂路径

osgeo.gdal.SetConfigOption('SHAPE_ENCODING','gb2312')#解决SHAPE⽂件的属性值

filename ='outlier_point.shp'

driver = ogr.GetDriverByName('ESRI Shapefile')

if os.access(filename,os.F_OK):

driver.DeleteDataSource(filename)

#以上，如果⽂件存在，就删除，重新创建

ds = driver.CreateDataSource(filename)

spatialref = osr.SpatialReference(proj_2)#空间参考

#如果直接写spatialref = proj_2会提⽰错误

#TypeError: in method 'DataSource_CreateLayer', argument 3 of type 'OSRSpatialReferenceShadow *' #因此要通过osr.SpatialReference将proj_2字符串转为类

#attention: 如果原空间参考是投影坐标那就是4666548,535646这种，如果

#原空间参考是经纬度坐标，则是12.4366655349731,-70.0561218261719这种

#不同的空间参考给出的原点处的位置参数也是不⼀样的

geomtype = ogr.wkbPoint #shape类型是点

layer = ds.CreateLayer(filename[:-4],srs = spatialref,geom_type = geomtype)

#layer = ds.CreateLayer(filename[:-4],srs = spatialref,geom_type = geomtype)

#layer图层的名称即.shp之前的字符串，空间参考是spatialref，⼏何对象是点

loc_p1_line = positive_outlier_index[0][0]#第⼀个点的⾏号

loc_p1_sample = positive_outlier_index[1][0]#列号

geo_p1_x = geotans_2[0]+loc_p1_sample*geotans_2[1]+ loc_p1_line*geotans_2[2]

geo_p1_y = geotans_2[3]+ loc_p1_sample*geotans_2[4]+ loc_p1_line*geotans_2[5]

pntp1_wkt ='POINT('+str(geo_p1_x)+' '+str(geo_p1_y)+')'

#attention：中间的是空格，不能是逗号

loc_p2_line = positive_outlier_index[0][1]

loc_p2_sample = positive_outlier_index[1][1]

geo_p2_x = geotans_2[0]+loc_p2_sample*geotans_2[1]+ loc_p2_line*geotans_2[2]

geo_p2_y = geotans_2[3]+ loc_p2_sample*geotans_2[4]+ loc_p2_line*geotans_2[5]

pntp2_wkt ='POINT('+str(geo_p2_x)+' '+str(geo_p2_y)+')'

loc_p3_line = positive_outlier_index[0][2]

loc_p3_sample = positive_outlier_index[1][2]

geo_p3_x = geotans_2[0]+ loc_p3_sample*geotans_2[1]+ loc_p3_line*geotans_2[2]

geo_p3_y = geotans_2[3]+ loc_p3_sample*geotans_2[4]+ loc_p3_line*geotans_2[5]

pntp3_wkt ='POINT('+str(geo_p3_x)+' '+str(geo_p3_y)+')'

loc_p4_line = positive_outlier_index[0][3]

loc_p4_sample = positive_outlier_index[1][3]

geo_p4_x = geotans_2[0]+ loc_p4_sample*geotans_2[1]+ loc_p4_line*geotans_2[2]

geo_p4_y = geotans_2[3]+ loc_p4_sample*geotans_2[4]+ loc_p4_line*geotans_2[5]

pntp4_wkt ='POINT('+str(geo_p4_x)+' '+str(geo_p4_y)+')'

loc_p5_line = positive_outlier_index[0][4]

loc_p5_sample = positive_outlier_index[1][4]

geo_p5_x = geotans_2[0]+ loc_p5_sample*geotans_2[1]+ loc_p5_line*geotans_2[2]

geo_p5_y = geotans_2[3]+ loc_p5_sample*geotans_2[4]+ loc_p5_line*geotans_2[5]

pntp5_wkt ='POINT('+str(geo_p5_x)+' '+str(geo_p5_y)+')'

理念英语#负异常值坐标求解

loc_n1_line = negative_outlier_index[0][0]

loc_n1_sample = negative_outlier_index[1][0]

geo_n1_x = geotans_2[0]+ loc_n1_sample*geotans_2[1]+ loc_n1_line*geotans_2[2]

geo_n1_y = geotans_2[3]+ loc_n1_sample*geotans_2[4]+ loc_n1_line*geotans_2[5]

pntn1_wkt ='POINT('+str(geo_n1_x)+' '+str(geo_n1_y)+')'

#-------------------------------write geomlist--------------------------

#geomlist是⼀个⾥⾯，⾥⾯存放着字符串，因此需要将此六个点放⼊这个list中

#geomlist存放的是⼏何形状，点线⾯的⼏何构成坐标

geomlist =[pntp1_wkt,pntp2_wkt,pntp3_wkt,pntp4_wkt,pntp5_wkt,pntn1_wkt]

#---------------------write field list-------------------------------

#field list 是个列表，列表⾥⾯是每⼀个feature的属性，是个字典型，如下：

#{'name': 'FIPS_CNTRY', 'type': 4, 'width': 2, 'decimal': 0}

#因此只需将异常值点字段信息按照列表中字典的形式写⼊这个列表供后⾯使⽤即可fieldlist =[{'name':'outlier_pt','type':4,'width':4,'decimal':0},\

{'name':'code','type':4,'width':2,'decimal':0},\

{'name':'color','type':4,'width':4,'decimal':0}]#创建字段list

#字段列表⾥每⼀个{}的内容即为shape属性表头⼀⾏的名称及其详细规定

#type: 0 int ；1 longinteger ；2 float ； 3 double； 4 text；5 date

#name 字段名称长度不能超过11个字母

#------------------------------------------------------------------

#本实验添加了三个字段，⼀个是“异常点”，⼀个是编码，⼀个是颜⾊

#⽤pi表⽰正异常值,i from 1 to n，代号1 ，颜⾊red

#⽤nj表⽰负异常值,j from 1 to m，代号-1，颜⾊blue

#python索引也是先⾏后列，⼤都是先⾏后列，只有IDL先

#---------------------write fieldlist into the lyer-------------------

for fd in fieldlist:

field = ogr.FieldDefn(fd['name'],fd['type'])

if'width'in fd:

field.SetWidth(fd['width'])

if'decimal'in fd:

field.SetPrecision(fd['decimal'])

layer.CreateField(field)

#-------------------------------write reclist-------------------------

#reclist为每个⼏何图形对应三个字段fieldlist的具体值

#如第⼀个图形outlierpoints叫p1，code为1，颜⾊是红⾊，是具体值

#如最后⼀个图形outlierpoints叫n1，code为-1，颜⾊是蓝⾊，是具体值

#本实验中的图形均为点

reclist =[{'outlier_pt':'p1','code':1,'color':'red'},\

{'outlier_pt':'p2','code':1,'color':'red'},\

{'outlier_pt':'p3','code':1,'color':'red'},\

{'outlier_pt':'p4','code':1,'color':'red'},\

{'outlier_pt':'p5','code':1,'color':'red'},\

{'outlier_pt':'n1','code':-1,'color':'blue'}]

#给出所创建好的图形(geomlist)中的字段(fieldlist)的具体值(reclist)

#----------------------------------write shape------------------------

for j in range(len(reclist)):

feat = ogr.Feature(layer.GetLayerDefn())

geom = ogr.CreateGeometryFromWkt(geomlist[j])

feat.SetGeometry(geom)

for f_d in fieldlist:

feat.SetField(f_d['name'],reclist[j][f_d['name']])

layer.CreateFeature(feat)

ds.Destroy()

为简化过程，将其中的geomlist写⼊过程和reclist写⼊过程写为循环，如下：（笔者很讨厌没有注释的代码，为保证代码可读性，依然保留了全部注释）

# -*- coding: utf-8 -*-失去的近义词

"""