Python空间数据处理环境搭建
 
作者：阿振
 
邮箱：tanzhenyugis@163.com
 
博客：https://blog.csdn.net/theonegis/article/details/80089375
 
修改时间：2018-04-26
 
优酷视频地址：http://v.youku.com/v_show/id_XMzU2NjAzMTk0NA==.html
 
 
Conda的下载和安装
 
什么是Conda? 官方定义：Package, dependency and environment management for any language—Python, R, Ruby, Lua, Scala, Java, JavaScript, C/ C++, FORTRAN
 
Conda就是一个虚拟环境和包（库）依赖管理工具
 
下载地址：Downloading conda
 
对于Windows版本的，确定Python版本和系统类型直接下载安装包进行安装
 
对于Linux和macOS系统，在Terminal中运行bash脚本进行安装即可。
 
Conda的使用
 
新建虚拟环境（Managing environments）
 
conda create -n osgeo
 
切换到新建的虚拟环境
 
source activate osgeo （Linux和macOS）
 
activate osgeo （Windows）
 
退出虚拟环境
 
source deactivate （Linux和macOS）
 
deactivate （Windows）
 
实用命令：
 
查看虚拟环境列表 conda env list 或者 conda info --envs
删除虚拟环境 conda remove --name <environment> --all
查看虚拟环境中的包列表 conda list
更新conda或者某个包 conda update conda 或者 conda update <package>
更新虚拟环境下的所有包 conda update --all
查看过时的包 conda search --outdated
搜索指定包 conda search <package>
删除某个包 conda remove <package>
添加channel到conda配置文件 conda config --add channels <channel> 或者 conda config --append channels <channel>
 
空间数据处理Python库的安装
 
常用的空间数据处理Python库
 
GDAL 全能型的基础空间数据处理库
fiona 基于GDAL的空间矢量数据处理库
rasterio 基于GDAL的空间栅格处理库
basemap 基于matplotlib的空间制图库
GeoPandas 基于pandas的空间数据分析库
RSGISlib 针对遥感数据及GIS分析的高级库
 
使用conda进行库的安装
 
打开命令行工具（Terminal），输入命令，进入虚拟环境
 
安装GDAL库
 
conda install -c conda-forge gdal
安装fiona库
 
conda install -c conda-forge fiona
安装rasterio库
 
conda install -c conda-forge rasterio
 
​
 
使用pip进行库的安装
 
什么是pip呢？pip是Python默认和推荐实用的包管理工具，可以利用pip从PyPI网络仓库自动下载Python包进行安装和管理。
 
对于Windows下的二进制库的预编译包，提供给大家一个网站：Unofficial Windows Binaries for Python Extension Packages
 
使用pip安装的时候，如果该Python包底层依赖一些C++库，则需要手动进行编译，或者安装指定平台下预编译好的库。
 
安装GDAL库
 
pip install GDAL‑2.2.4‑cp37‑cp37m‑win_amd64.whl
安装fiona库
 
pip install Fiona‑1.7.11.post1‑cp37‑cp37m‑win_amd64.whl
安装rasterio库
 
pip install rasterio‑1.0a12‑cp37‑cp37m‑win_amd64.whl
安装Jupyter
 
pip install jupyter
 
使用Jupyter Notebook进行编程
 
jupyter notebook

GDAL是空间数据处理的开源包，支持多种数据格式的读写。遥感图像是一种带大地坐标的栅格数据，遥感图像的栅格模型包含以下两部分的内容：
 
栅格矩阵：由正方形或者矩形栅格点组成，每个栅格点所对应的数值为该点的像元值，在遥感图像中用于表示地物属性值；遥感图像有单波段与多波段，波段表示地物属性的种类，每个波段表示地物一种属性。
 
大地坐标：空间数据参考表示地图的投影信息；仿射矩阵能将行列坐标映射到面坐标上。
 
GDAL读写遥感数据的代码：
 
from osgeo import gdal
import os

class GRID:

    #读图像文件
    def read_img(self,filename):
        dataset=gdal.Open(filename)       #打开文件

        im_width = dataset.RasterXSize    #栅格矩阵的列数
        im_height = dataset.RasterYSize   #栅格矩阵的行数

        im_geotrans = dataset.GetGeoTransform()  #仿射矩阵
        im_proj = dataset.GetProjection() #地图投影信息
        im_data = dataset.ReadAsArray(0,0,im_width,im_height) #将数据写成数组，对应栅格矩阵

        del dataset 
        return im_proj,im_geotrans,im_data

    #写文件，以写成tif为例
    def write_img(self,filename,im_proj,im_geotrans,im_data):
        #gdal数据类型包括
        #gdal.GDT_Byte, 
        #gdal .GDT_UInt16, gdal.GDT_Int16, gdal.GDT_UInt32, gdal.GDT_Int32,
        #gdal.GDT_Float32, gdal.GDT_Float64

        #判断栅格数据的数据类型
        if 'int8' in im_data.dtype.name:
            datatype = gdal.GDT_Byte
        elif 'int16' in im_data.dtype.name:
            datatype = gdal.GDT_UInt16
        else:
            datatype = gdal.GDT_Float32

        #判读数组维数
        if len(im_data.shape) == 3:
            im_bands, im_height, im_width = im_data.shape
        else:
            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

if __name__ == "__main__":
    os.chdir(r'D:\Python_Practice')                        #切换路径到待处理图像所在文件夹
    run = GRID()
    proj,geotrans,data = run.read_img('LC81230402013164LGN00.tif')        #读数据
    print proj
    print geotrans
    print data
    print data.shape
    run.write_img('LC81230402013164LGN00_Rewrite.tif',proj,geotrans,data) #写数据
 
在GDAL遥感影像读写的基础上，我们可以进行遥感图像的各种公式计算和统计分析。 
 例如我们所熟知的计算NDVI（归一化植被指数），只要在以上代码倒数第二行中插入代码：
 
import numpy as np
data = data.astype(np.float)
ndvi = (data[3]-data[2])/(data[3]+data[2])                         #3为近红外波段；2为红波段
run.write_img('LC81230402013164LGN00_ndvi.tif',proj,geotrans,ndvi) #写为ndvi图像
 
当然，这是理想的NDVI，实际处理NDVI还会遇到一些其他要处理的问题。例如NDVI值应该在区间[-1，1]内，但实际中会出现大于1或小于-1的情况，或者某些像点是坏点，出现空值nan，需要进一步的配套处理。

HDF格式 
 一种能有条不紊 、完备地保存遥感影像的属性和空间信息数据 , 同时使查询和提取相关数据也很方便容易的数据格式格式。HDF格式是一种较为常见的遥感数据格式（例如MODIS、3B43降水数据），但GDAL暂时未能给HDF格式提供支持，需要PYHDF与GDAL交互，从而将HDF格式转换为GDAL支持的格式。
 
PyHDF 
 python处理HDF数据的开源包，笔者安装的pyhdf4(电脑配置WIN64+python2.7)，在我的另一篇文章《Python开源包安装：PYHDF安装》中，有详细介绍PyHDF的安装方法。
 
下面笔者就以 MODIS 1B 数据为例，根据PyHDF和GDAL，介绍转HDF格式影像为包含红、红外两个波段的GDAL所支持的栅格格式影像。
 
1、pyhdf 获取待处理影像栅格矩阵 
 导入需要用到的模块：
 
import os
from pyhdf.SD import SD
import numpy as np
 
将工作空间转到待处理影像所在文件夹：
 
os.chdir(r'D:\modis_ndvi_classificate\fix_proj\pre_ndvi')
 
SD打开待处理影像
 
FILE_NAME = 'MOD09Q1.A2015169.h26v04.006.2015301111712.hdf'
hdf = SD(FILE_NAME, SDC.READ)
 
查看影像内所包含的数据：
 
print hdf.datasets()
 
输入“print hdf.datasets()”后，显示： 
  
 其中，sur_refl_b02, sur_refl_b01 分别代表MODIS的红外、红两个波段。
 
获取红、红外两个波段的影像栅格，data1与data2：
 
data1 = hdf.select('sur_refl_b01')
data1 = data1[:]
data2 = hdf.select('sur_refl_b02')
data2 = data2[:]
 
将data1与data2叠加：
 
data = np.concatenate ([data1[np.newaxis,:], data2[np.newaxis,:]], axis=0)
 
2、获取仿射矩阵 
 查看影像头文件信息：
 
attr = hdf.attributes(full=1)
attNames = attr.keys()
attNames.sort()
print attNames
 
输入”print attNames” 后会显示该影像所包含头文件的目录有[‘ArchiveMetadata.0’, ‘CoreMetadata.0’,’HDFEOSVersion’,’StructMetadata.0’,’identifier_product_doi’,’identifier_product_doi_authority’]，我们打开’ArchiveMetadata.0’，中间包含了影像的起始点坐标和分辨率信息，：
 
t = attr['ArchiveMetadata.0']
print t[0]
 
输入‘print t[0]’后，我们获得该影像起始点坐标：经度104.43258313171，纬度49.9999999955098；红，红外波段分辨率为231.656358263889，则该影像仿射矩阵为：
 
geot = (104.43258313171,231.656358263889,0,49.9999999955098,0,-231.656358263889)
 
3、给影像添加用户所指定的投影： 
 （以WGS84为例）
 
from osgeo import osr
proj= osr.SpatialReference()
proj.ImportFromEPSG(4326)
proj.ExportToWkt()
 
4、重写影像 
 有了影像的栅格矩阵、仿射矩阵，投影信息后，我们就可以按照《Python空间数据处理1： GDAL读写遥感图像》部分的内容，重新将影像写成GDAL所支持的栅格格式文件了。

 
流程图如下：
 
 
 
用到的函数：
 （1）object=ogr.SpatialReference()
 该函数是SpatialReference类的一个构造函数，它创建了一个SpatialReference类的对象，用以存储空间参考信息。
 （2）object.ImportFromEPSG(4326)
 该函数是SpatialReference类的一个成员函数，它从EPSG中导入一个空间参考。object.ImportFromEPSG(4326)表示从EPSG中导入WGS84参考（点击此处可以查看EPSG里面各种空间参考以及编号）。
 （3） trans=osr.CoordinateTransformation(spatialref_source,spatialref_target)
 该函数是CoordinateTransformation类（该类用于在不同坐标系之间进行转换）的构造函数，构建一个CoordinateTransformation对象，接下来可以利用该对象访问CoordinateTransformation类的成员函数以完成点坐标在不同空间参考下的转换。
 （4）after_trans_point=trans.TransformPoint(x,y)
 该函数是CoordinateTransformation类的一个成员函数，使用CoordinateTransformation类的对象调用该函数可以完成点的坐标在不同空间参考下的转换。该函数的返回值为元祖（tuple）。
 （5）after_trans_points=trans.TransformPoints([(x1,y1),(x2,y2)])
 该函数同样为CoordinateTransformation类的成员函数，用来转换多个点，该函数的形式参数是一个列表。
 
 
点的转换代码
 
 
#coding=utf-8 
from osgeo import osr,ogr,gdal
def point_transform(source_ref,target_ref,x,y):
#创建目标空间参考
    spatialref_target=osr.SpatialReference()
    spatialref_target.ImportFromEPSG(target_ref) #西安80 高斯可吕格投影
    #创建原始空间参考
    spatialref_source=osr.SpatialReference()
    spatialref_source.ImportFromEPSG(source_ref)  #WGS84
        #构建坐标转换对象，用以转换不同空间参考下的坐标
    trans=osr.CoordinateTransformation(spatialref_source,spatialref_target)
    coordinate_after_trans=trans.TransformPoint(x,y)
    return coordinate_after_trans
	 #以下为转换多个点（要使用list）
    #coordinate_trans_points=trans.TransformPoints([(117,40),(120,36)])
    #print(coordinate_trans_points)


#函数调用
point_trans=point_transform(4326,2331,120,36)
#4326 为原始空间参考的ESPG编号  2331为目标空间参考的ESPG编号
#120为经度，36为纬度
print(point_trans)

 
注：#coding=utf-8 用来解决代码中有中文时会出错的问题
 
结果