python gdal 基于栅格shp文件裁剪geotif图

import numpy as np
import gdal
import operator
from functools import reduce
import os
import shapefile

def geo2imagexy(dataset, x, y):
‘’’
根据GDAL的六参数模型将给定的投影或地理坐标转为影像图上坐标
:param dataset: GDAL地理数据
:param x: 投影或地理坐标x
:param y: 投影或地理坐标y
:return: 影坐标或地理坐标(x, y)对应的影像图上坐标
‘’’
trans = dataset.GetGeoTransform()
a = np.array([[trans[1], trans[2]], [trans[4], trans[5]]])
b = np.array([x - trans[0], y - trans[3]])
return np.linalg.solve(a, b) # 使用numpy的linalg.solve进行二元一次方程的求解

#读取要裁剪的原图
in_ds = gdal.Open(“raw.tif”)
print(“open tif file succeed”)

#读取原图中的每个波段
in_band1 = in_ds.GetRasterBand(1)
in_band2 = in_ds.GetRasterBand(2)
in_band3 = in_ds.GetRasterBand(3)
in_band4 = in_ds.GetRasterBand(4)

#导入栅格shp文件
shp_path=‘D:/’
r = shapefile.Reader(shp_path+‘1.shp’)
minX, minY, maxX, maxY = r.bbox

#定义切图的起始点和终点坐标(相比原点的横坐标和纵坐标偏移量)
offset_x,offset_y=geo2imagexy(in_ds,minX,minY)
endset_x,endset_y=geo2imagexy(in_ds,maxX,maxY)

#定义切图的大小（矩形框）
block_xsize=int(endset_x-offset_x+1)
block_ysize=int(abs(endset_y-offset_y)+1)

#从每个波段中裁剪需要的矩形框内的数据
out_band1 = in_band1.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
out_band2 = in_band2.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
out_band3 = in_band3.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)
out_band4 = in_band4.ReadAsArray(offset_x, offset_y, block_xsize, block_ysize)

#获取Tif的驱动，为创建切出来的图文件做准备
gtif_driver = gdal.GetDriverByName(“GTiff”)

#创建切出来的要存的文件（4代表4个不都按，最后一个参数为数据类型，跟原文件一致）
out_ds = gtif_driver.Create(‘clip.tif’, block_xsize, block_ysize, 4, in_band1.DataType)
print(“create new tif file succeed”)

#获取原图的原点坐标信息
ori_transform = in_ds.GetGeoTransform()

#读取原图仿射变换参数值
top_left_x = ori_transform[0] # 左上角x坐标
w_e_pixel_resolution = ori_transform[1] # 东西方向像素分辨率
top_left_y = ori_transform[3] # 左上角y坐标
n_s_pixel_resolution = ori_transform[5] # 南北方向像素分辨率

#根据反射变换参数计算新图的原点坐标
top_left_x = top_left_x + offset_x * w_e_pixel_resolution
top_left_y = top_left_y + offset_y * n_s_pixel_resolution

#将计算后的值组装为一个元组，以方便设置
dst_transform = (top_left_x, ori_transform[1], ori_transform[2], top_left_y, ori_transform[4], ori_transform[5])

#设置裁剪出来图的原点坐标
out_ds.SetGeoTransform(dst_transform)

#设置SRS属性（投影信息）
out_ds.SetProjection(in_ds.GetProjection())

#写入目标文件
out_ds.GetRasterBand(1).WriteArray(out_band1)
out_ds.GetRasterBand(2).WriteArray(out_band2)
out_ds.GetRasterBand(3).WriteArray(out_band3)
out_ds.GetRasterBand(4).WriteArray(out_band4)

#将缓存写入磁盘
out_ds.FlushCache()
print(“FlushCache succeed”)
del out_ds
print(“End!”)

python gdal 基于栅格shp文件裁剪geotif图相关推荐

python shp文件_对python 读取线的shp文件实例详解
如下所示: import shapefile sf = shapefile.reader("e:\\1.2\\cs\\dx_csl.shp") shapes = sf.shapes ...
Python 画降雨图，mash 数据掩膜裁剪出想要的区域，根据shp文件裁剪tif数据
1. 背景介绍 2. 导库 import numpy as np import xarray as xr import matplotlib.pyplot as plt from matplotlib ...
NC文件不规则裁剪(利用shp文件裁剪)
文章目录示例数据与软件介绍数据软件中国区域提取利用Python裁剪nc文件示例数据与软件介绍数据 NCEP的逐月2m气温数据中国基础shp文件软件这次除了常用的Python之外,我 ...
gdal java shp_【GDAL/OGR】利用GDAL/OGR读取shp文件并转换为json文件（Java版）
前言: 对于GIS开发者来说,GDAL/OGR是最熟悉不过的开源GIS库了,GDAL(Geospatial Data Abstraction Library)是一个在X/MIT许可协议下的开源栅格空间 ...
Python gdal库读取tif文件
from osgeo import gdal # GDAL库主要提供对栅格数据的处理,使用抽象数据模型来解析所支持的数据格式 import filename_cut as fc import matp ...
python文件的两种类型进制文件,Python之基于十六进制判断文件类型
#!/usr/bin/env python#-*- coding: utf-8 -*-#@Author : suk importstructfrom io importBytesIO#支持文件类型#用 ...
Python中ArcPy栅格裁剪栅格：批量对齐栅格图像范围并统一行数与列数
本文介绍基于Python中ArcPy模块,实现基于栅格图像批量裁剪栅格图像,同时对齐各个栅格图像的空间范围,统一其各自行数与列数的方法. 首先明确一下我们的需求.现有某一地区的多张栅格遥感影像 ...
Android GIS开发系列-- 入门季（13）Gdal简单写个shp文件
Gdal是用来读写栅格与矢量数据的,在Gdal官网,可以下载相关的资源进行平台的编译.其实Arcgis底层也是用Gdal来读取shp文件的,那在Android中可以直接读写shp文件吗,是可以的.这里 ...
Arcgis加载shp文件
Arcgis地图加载shp文件,效果如图(当前shp文件包含雨量等值面数据): 1.思路:使用js-shapefile-to-geojson把shp和dbf文件转为geojson格式数据,然后通过Ar ...
python使用gdal将shp文件转为TIF
python使用gdal将shp文件转为TIF 方法一 # 缺少获取shp文件坐标系的步骤 def vector2raster(inputfilePath, outputfile, resp):sf ...

python gdal 基于栅格shp文件裁剪geotif图

python gdal 基于栅格shp文件裁剪geotif图相关推荐

最新文章

热门文章