python处理shp和栅格文件的相关库shapefile、gdal等

读取含polygon的shp文件:

def readshp(shp_path):
	sf = shapefile.Reader(shp_path)#创建reader类的对象进行shapefile文件的读取
	shapes = sf.shapes()# .shapes()读取几何数据信息，存放着该文件中所有对象的 几何数据
	#records = sf.records()
	out = []
	for i, shape in enumerate(Shapes):
		polybox = Polygon(shape.points)#shape.points=[(x1,y1),(x2,y2),(x3,y3),…],首尾相连，起始点重复，矩形五个点坐标处于一个列表
		out.append(polybox)
	return out

下图是该函数调试一张含多个矩形polygon的shp文件结果

def read_TIF(filename):
	dataset = gdal.Open(filename)
	img_width = dataset.RasterXSize#获取栅格宽度 水平方向
	img_height = dataset.RasterYSize#获取栅格高度 竖直方向
	img_nbands = dataset.RasterCount#获取栅格通道数
	geoTransform = dataset.GetGeoTransform()
	x0, y0 = geoTransform[0], geoTransform[3]
	x1, y1 = geoTransform[0] + img_width * geoTransform[1], geoTransform[3] + img_height * geoTransform[5]
	polybox_tif = Polygon([(x0,y0),(x1,y0),(x1,y1),(x0,y1)])
	image_data = dataset.ReadAsArray(0,0,image_width,image_height)
	img_array = img_data.reshape((img_height,img_width,4))
	......

其中，对于gdal库，GDAL中使用dataset表示一个栅格数据(使用抽象类GDALDataset表示)，一个dataset包含了对于栅格数据的波段，空间参考以及元数据等信息。一张GeoTIFF遥感影像，一张DEM影像，或者一张土地利用图，在GDAL中都是一个GDALDataset。

class Grid:
    # 读取图像文件
    def read_img(self, filename):
        data = gdal.Open(filename)  # 打开文件
        im_width = data.RasterXSize  # 读取图像行数
        im_height = data.RasterYSize  # 读取图像列数
        im_geotrans = data.GetGeoTransform()
        # 仿射矩阵，左上角像素的大地坐标和像素分辨率。
        # 共有六个参数，分表代表左上角x坐标；东西方向上图像的分辨率；如果北边朝上，地图的旋转角度，0表示图像的行与x轴平行；左上角y坐标；
        # 如果北边朝上，地图的旋转角度，0表示图像的列与y轴平行；南北方向上地图的分辨率。
        im_proj = data.GetProjection()  # 地图投影信息
        im_data = data.ReadAsArray(0, 0, im_width, im_height)  # 此处读取整张图像
        # ReadAsArray(, , , )
        # 读出从(xoff,yoff)开始，大小为(xsize,ysize)的矩阵。
        del data  # 释放内存，如果不释放，在arcgis，envi中打开该图像时会显示文件被占用
        return im_proj, im_geotrans, im_data

对于GetGeoTransform()
GetGeoTransform是仿射变换，计算坐标偏移。（现实世界坐标（地理坐标或投影坐标）与图像坐标行列号），其中Geo[0/3]是现实世界坐标而不是行列号。
0、3 是x和y坐标 起始点现实世界坐标 1、5 是像素宽度和高度 2、4 是x和y方向旋转角

#GetGeoTransform，在真实坐标和栅格数据坐标具有相同srs情况下，计算坐标偏移。
#作用：图像坐标（行列号）和现实世界坐标（投影坐标或地理坐标）之间的转换。是仿射变换，不是投影转换，和上面不同。
#0、3 x\y坐标 起始点现实世界坐标  1、5 像素宽度和高度  2、4 x\y方向旋转角
gt = ds.GetGeoTransform() #正变换：图像坐标到现实世界坐标。正变换时输入行列号，输出的现实世界坐标是栅格图像srs下的坐标
inv_gt = gdal.InvGeoTransform(gt) #逆变换：现实世界坐标到图像坐标
offsets = gdal.ApplyGeoTransform(inv_gt, 465200, 5296000) #逆变换：输入的投影坐标具有和栅格图像的相同的srs
xoff, yoff = map(int, offsets)  #取整