基于模板测试实现半透明多边形运算 .

多边形常见操作

android api sample中的region 操作如下图:


本文基于opengl 的stencil buffer 实现这几种多边形运算。

OpenGL 模板测试实现

如下图:

 如上图,
Source为原始半透明重叠多边形绘制效果,重叠部分因为融合的缘故 有增强效果。

Union   联合操作,需保证每个像素 当且仅当只属于一个多边形,为此通过模板测试实现的具体思路:
  1. 用0x00 清除模板缓冲区,开启模板测试
  2. 写模板缓冲区,同时进行模板测试,同时写颜色缓冲区:保证只有模板值为0的片源可以通过模板测试,而且通过模板测试后 模板值加1 操作。

Difference 做差操作,是两个矩形相减效果,顺序无关,具体实现:
  1. 0x00 清除模板缓冲区
  2. 禁止颜色缓冲区写入,禁止模板测试,画任意一个矩形,将模板值加1 操作。此时末班缓冲区中第一个矩形绘制区域为mask区域,该区域内不能绘图。
  3. 开启模板测试,只能在模板值为0x00的地方绘图,画第二个矩形,模板测试通过后模板值加1。

XOr 异或操作,两个矩形相交的区域不显示,不相交的区域显示,具体思路:
  1. 0x00 清除模板缓冲区,
  2. 禁掉模板测试,禁掉颜色缓冲区写入,将模板值加1操作,绘制两个矩形,相交区域的模板值为2,其他矩形区域模板值为1。
  3. 开启模板测试,允许写入颜色缓冲区,测试条件为模板值等于1,模板值为8位,此时掩码为0x01,二进制为00000001。绘制两个矩形。模板测试通过后将模板值置为0。

Intersection 相交操作,只显示两个矩形相交区域,具体思路:
  1. 0x00 清除模板缓冲区,
  2. 禁掉模板测试,禁掉颜色缓冲区写入,将模板值加1操作,绘制两个矩形,相交区域的模板值为2,其他矩形区域模板值为1。
  3. 开启模板测试,允许写入颜色缓冲区,测试条件为模板值等于2,模板值为8位,此时掩码为0x03,二进制为00000011。绘制任意一个矩形。模板测试通过后将模板值置为0。

GLUT框架下实现:

// @module:		stencil render
// @file:		opacityRectRender.cpp
// @author:		peteryfren
// @date:		2013/3/28
// @version:	1.0
// @desc:		通过模板缓冲区实现半透明矩形的常见操作。

#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <assert.h>
#include <time.h>
#include <math.h>

#include <gl/glut.h>

const int windowWidth = 256;
const int windowHeight = 256;

GLenum checkForError(char *loc);

void opacitySourceRender(int width, int height);
void opacityRectXorRender(int width, int height);
void opacityRectUnionRender(int width, int height);
void opacityRectIntersectRender(int width, int height);
void opacityRectDifferenceRender(int width, int height);

void Redraw(void)
{
	glCullFace(GL_BACK);

	float destAlpha = 0.5;
	glClearColor(0,0,0, 0.5);	// 指定清理颜色值
	glClearStencil(0x00);		// 指定清理帧缓冲值
	glClearDepth(0.0);			// 指定清理深度缓冲值
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT | GL_STENCIL_BUFFER_BIT);

	// glEnable (GL_TEXTURE_2D);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

	opacitySourceRender(windowWidth, windowHeight);
	
	// opacityRectDifferenceRender(windowWidth, windowHeight);
	// opacityRectIntersectRender(windowWidth, windowHeight);
	// opacityRectUnionRender(windowWidth, windowHeight);
	// opacityRectXorRender(windowWidth, windowHeight);

	glutSwapBuffers();

	checkForError("swap");
}

void opacitySourceRender(int width, int height)
{
	// one way to solve this in generic OpenGL is to 
	//		fill the framebuffer with your desired alpha channel, 
	// switch the blending mode to glBlendFunc(GL_ONE_MINUS_DST_ALPHA, GL_DST_ALPHA)
	//		and draw the rectangles. 

	glBlendFunc(GL_ONE_MINUS_DST_ALPHA, GL_DST_ALPHA);

	float tileVertices[] = {0,height/2, width,height/2, width,height*3/4, 0,height*3/4}; 
	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	float tileVertices1[] = {width/4,0, width/2,0, width/2,height, width/4,height}; 
	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices1);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);
}

void opacityRectDifferenceRender(int width, int height)
{
	glEnable(GL_STENCIL_TEST);

	// draw mask area
	glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE);
	glStencilFunc(GL_ALWAYS, 0, 0xFF);			// 禁模板测试
	glStencilOp(GL_KEEP, GL_KEEP, GL_INCR);		// 模板值加1

	glColor4f(1, 0, 0, 0.1);
	float tileVertices[] = {0,height/2, width,height/2, width,height*3/4, 0,height*3/4};
	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);

	//static unsigned char data[windowWidth*windowHeight];
	//glReadPixels(0, 0, windowWidth, windowHeight, GL_STENCIL_INDEX, GL_UNSIGNED_BYTE, data);
	//SaveFileGrayBMP("stencil.bmp", data, windowWidth,  windowHeight);

	// draw main area
	glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
	glStencilFunc(GL_EQUAL, 0x00, 0x01);	// 只能在模板为0的地方绘图
	glStencilOp(GL_KEEP, GL_KEEP, GL_INCR);	

	glColor4f(1, 0, 0, 0.5);
	float tileVertices1[] = {width/4,0, width/2,0, width/2,height, width/4,height}; 
	glVertexPointer(2, GL_FLOAT, 0, tileVertices1);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	glDisable(GL_STENCIL_TEST);
}

void opacityRectXorRender(int width, int height)
{
	glEnable(GL_STENCIL_TEST);

	glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE);

	// draw mask area 
	glStencilOp(GL_INCR, GL_INCR, GL_INCR);
	glStencilFunc(GL_ALWAYS, 0x00, 0xFF);		// 禁用模板

	float tileVertices[] = {0,height/2, width,height/2, width,height*3/4, 0,height*3/4};
	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);

	glColor4f(1, 0, 0, 0.5);
	float tileVertices1[] = {width/4,0, width/2,0, width/2,height, width/4,height}; 
	glVertexPointer(2, GL_FLOAT, 0, tileVertices1);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	//static unsigned char data[windowWidth*windowHeight];
	//glReadPixels(0, 0, windowWidth, windowHeight, GL_STENCIL_INDEX, GL_UNSIGNED_BYTE, data);
	//SaveFileGrayBMP("stencil.bmp", data, windowWidth,  windowHeight);

	glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
	glStencilFunc(GL_EQUAL, 0x01, 0x01);		// 模板值为1的区域才绘制
	glStencilOp(GL_ZERO, GL_ZERO, GL_ZERO);		// 用0替换当前值

	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices1);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	glDisable(GL_STENCIL_TEST);
}

void opacityRectIntersectRender(int width, int height)
{
	glEnable(GL_STENCIL_TEST);
	
	glColorMask(GL_FALSE, GL_FALSE, GL_FALSE, GL_FALSE);
	
	// draw mask area 
	glStencilOp(GL_INCR, GL_INCR, GL_INCR);
	glStencilFunc(GL_ALWAYS, 0x00, 0xFF);		// 禁用模板

	float tileVertices[] = {0,height/2, width,height/2, width,height*3/4, 0,height*3/4};
	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);

	glColor4f(1, 0, 0, 0.5);
	float tileVertices1[] = {width/4,0, width/2,0, width/2,height, width/4,height}; 
	glVertexPointer(2, GL_FLOAT, 0, tileVertices1);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	//static unsigned char data[windowWidth*windowHeight];
	//glReadPixels(0, 0, windowWidth, windowHeight, GL_STENCIL_INDEX, GL_UNSIGNED_BYTE, data);
	//SaveFileGrayBMP("stencil.bmp", data, windowWidth,  windowHeight);

	glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);
	glStencilFunc(GL_EQUAL, 0x02, 0x03);		// 模板值为2的区域才绘制 掩码0000 0011
	glStencilOp(GL_ZERO, GL_ZERO, GL_ZERO);

	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices1);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	glDisable(GL_STENCIL_TEST);
}

void opacityRectUnionRender(int width, int height)
{
	// 模板缓冲区就一个字节

	glEnable (GL_STENCIL_TEST);		

	// glStencilFunc(GL_ALWAYS,0,0xFF);	// 禁用模板测试

	glColorMask(GL_TRUE, GL_TRUE, GL_TRUE, GL_TRUE);	// 开启写入颜色缓冲区

	// fail = KEEP  模版测试失败后保持模板缓冲区中值 不变
	// zfail =  KEEP 深度测试失败后保持模版缓冲区中值 不变
	// zpass = INCR 深度测试通过后对模版缓冲区中值 加1操作
	glStencilOp(GL_KEEP, GL_KEEP, GL_INCR);		
	
	// glStencilOp(GL_KEEP, GL_KEEP, GL_KEEP);

	// func = EQUAL 只有与参考值相等才能通过模版测试
	// ref = 0x00	参考值
	// mask = 0x01	只比较最低位
	glStencilFunc(GL_EQUAL, 0x0, 0x1);	// 开启模板测试

	glEnable (GL_BLEND);		// 开启混合
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);

	// sfactor = 1 - dstA
	// dfactor = dstA
	// glBlendFunc(GL_ONE_MINUS_DST_ALPHA, GL_DST_ALPHA);

	float tileVertices[] = {0,height/2, width,height/2, width,height*3/4, 0,height*3/4};
	glColor4f(1, 0, 0, 0.5);
	glVertexPointer(2, GL_FLOAT, 0, tileVertices);
	glEnableClientState(GL_VERTEX_ARRAY);
	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);

	glColor4f(1, 0, 0, 0.5);
	float tileVertices1[] = {width/4,0, width/2,0, width/2,height, width/4,height}; 
 	glVertexPointer(2, GL_FLOAT, 0, tileVertices1);
 	glEnableClientState(GL_VERTEX_ARRAY);
 	glDrawArrays(GL_TRIANGLE_FAN, 0, 4);
	glColor4f(1, 1, 1, 1);

	//static unsigned char data[windowWidth*windowHeight];
	//glReadPixels(0, 0, windowWidth, windowHeight, GL_STENCIL_INDEX, GL_UNSIGNED_BYTE, data);
	//SaveFileGrayBMP("stencil.bmp", data, windowWidth,  windowHeight);

	glDisable (GL_STENCIL_TEST);
}

void Button(int button, int state, int x, int y)
{
	if (state != GLUT_UP)
		return;

	switch (button) 
	{
	case GLUT_LEFT_BUTTON:

		glutPostRedisplay();
		break;
	
	case GLUT_RIGHT_BUTTON:

		glutPostRedisplay();
		break;
	}
}

void Keyboard(unsigned char key, int x, int y)
{

	glutPostRedisplay();
}

void Reshape(int width, int height)
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(0.0f, (GLfloat) windowWidth, 0.0f, 
		(GLfloat) windowHeight, -1.0f, 1.0f);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glViewport(0, 0, windowWidth, windowHeight);
}

int main(int argc, char *argv[])
{
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_RGB|GLUT_DOUBLE|GLUT_DEPTH|GLUT_STENCIL);

	glutInitWindowSize(windowWidth, windowHeight);
	glutCreateWindow("Region");

	// set up world space to screen mapping
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(0.0f, (GLfloat) windowWidth, 0.0f, 
		(GLfloat) windowHeight, -1.0f, 1.0f);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glViewport(0, 0, windowWidth, windowHeight);

	glClearColor(1.0, 1.0, 1.0, 1.0);
	glutDisplayFunc(Redraw);
	glutReshapeFunc(Reshape);
	glutMouseFunc(Button);
	glutKeyboardFunc(Keyboard);

	glEnable(GL_CULL_FACE);

	glutMainLoop();

	return 0;
}

GLenum checkForError(char *loc)
{
	GLenum errCode;
	const GLubyte *errString;

	if ((errCode = glGetError()) != GL_NO_ERROR)
	{
		errString = gluErrorString(errCode);
		printf("OpenGL error: %s",errString);

		if (loc != NULL)
			printf("(%s)",loc);

		printf("\n");
	}

	return errCode;
}


 

参考

3. opencsg   http://www.opencsg.org/

 

1. 程序功能概述 (先总体介绍程序具有哪些功能,再依次介绍每一项功能。) 1.1 功能xxx概述 1.2 功能xxx概述 2. 功能实现原理及方法 2.1 功能xxx实现原理及方法 有关键的程序设计思路与方法介绍 2.2 功能xxx实现原理及方法 有关键的程序设计思路与方法介绍 3. 结果展示与说明 3.1功能xxx演示结果与说明 有程序运行结果截图与说明 3.2 功能xxx演示结果与说明 有程序运行结果截图与说明 4. 心得体会 上面这个是我的大作业的程序说明文档的模板,请你根据我的程序以及代码写一个程序说明文档,然后我的功能包括打开矢量数据,打开栅格数据,矢量数据的分析功能包括缓冲区分析,属性查询以及绘制质心,空间查询,保存矢量数据,然后栅格数据的功能包括波段组合,栅格裁剪,以及波段运算,还有一个保存栅格数据,然后这个是我的代码# main_window.py(主窗口逻辑) import os import numpy as np from PySide6.QtWidgets import QMainWindow, QFileDialog, QGraphicsScene, QGraphicsView, QMessageBox, QGraphicsPathItem from PySide6.QtGui import QPainterPath, QPen, QBrush, QAction, QTransform, QImage, QPixmap, QColor from PySide6.QtCore import Qt, QRectF, QPointF from osgeo import ogr, gdal from PySide6.QtWidgets import QInputDialog # 新增输入对话框 # 新增自定义图形项类(用于存储属性) class FeatureItem(QGraphicsPathItem): def __init__(self, path, attributes): super().__init__(path) self.attributes = attributes # 存储属性字典 class MainWindow(QMainWindow): def __init__(self): super().__init__() self.setWindowTitle("GIS软件") self.setGeometry(100, 100, 800, 600) ogr.UseExceptions() self.init_ui() self.scene = QGraphicsScene(self) self.graphicsView.setScene(self.scene) # 新增:存储所有几何边界 self.total_bounds = QRectF() def init_ui(self): self.toolBar = self.addToolBar("工具") self.actionOpen_Vector_Data = QAction("打开矢量数据", self) self.toolBar.addAction(self.actionOpen_Vector_Data) # 新增栅格动作 self.actionOpen_Raster_Data = QAction("打开栅格数据", self) self.toolBar.addAction(self.actionOpen_Raster_Data) # 添加到工具栏 # 新增缓冲区分析按钮 self.actionBuffer_Analysis = QAction("缓冲区分析", self) self.toolBar.addAction(self.actionBuffer_Analysis) self.graphicsView = QGraphicsView() self.setCentralWidget(self.graphicsView) # 新增属性查询按钮 self.actionQuery_Attribute = QAction("属性查询", self) self.toolBar.addAction(self.actionQuery_Attribute) self.actionOpen_Vector_Data.triggered.connect(self.open_vector_data) self.actionOpen_Raster_Data.triggered.connect(self.open_raster_data) # 新增连接 self.actionBuffer_Analysis.triggered.connect(self.buffer_analysis) self.actionQuery_Attribute.triggered.connect(self.enable_query_mode) # 新增鼠标点击事件 self.graphicsView.setMouseTracking(True) self.is_query_mode = False # 新增波段组合按钮 self.actionBand_Combination = QAction("波段组合", self) self.toolBar.addAction(self.actionBand_Combination) self.actionBand_Combination.triggered.connect(self.open_band_combination) # 新增栅格裁剪按钮(在init_ui方法末尾添加) self.actionClip_Raster = QAction("栅格裁剪", self) self.toolBar.addAction(self.actionClip_Raster) self.actionClip_Raster.triggered.connect(self.clip_raster) # 新增连接 self.actionBand_Calculation = QAction("波段运算", self) self.toolBar.addAction(self.actionBand_Calculation) self.actionBand_Calculation.triggered.connect(self.band_calculation) # 新增质心绘制按钮(放在init_ui方法中) self.actionDraw_Centroids = QAction("绘制质心", self) self.toolBar.addAction(self.actionDraw_Centroids) self.actionDraw_Centroids.triggered.connect(self.draw_centroids) self.centroid_items = [] # 新增:存储质心图形项 # 新增空间查询按钮(放在init_ui方法中) self.actionSpatial_Query = QAction("空间查询", self) self.toolBar.addAction(self.actionSpatial_Query) self.actionSpatial_Query.triggered.connect(self.enable_spatial_query_mode) self.is_spatial_query_mode = False self.spatial_query_results = [] # 存储查询结果 # 新增保存按钮(放在init_ui方法中) self.actionSave_Vector = QAction("保存矢量数据", self) self.toolBar.addAction(self.actionSave_Vector) self.actionSave_Vector.triggered.connect(self.save_vector_data) # 新增栅格保存按钮(放在init_ui方法中) self.actionSave_Raster = QAction("保存栅格数据", self) self.toolBar.addAction(self.actionSave_Raster) self.actionSave_Raster.triggered.connect(self.save_raster_data) def open_vector_data(self): file_path, _ = QFileDialog.getOpenFileName( self, "打开矢量文件", "", "Shapefile (*.shp);;GeoJSON (*.geojson);;All Files (*)" ) if file_path: self.load_vector_data(file_path) # 新增:自动缩放视图 self.auto_zoom() def load_vector_data(self, file_path): self.scene.clear() self.total_bounds = QRectF() # 重置边界 try: data_source = ogr.Open(file_path, 0) layer = data_source.GetLayer(0) for feature in layer: geom = feature.GetGeometryRef() path = self.geometry_to_qpainterpath(geom) # 更新总边界 if path.boundingRect().isValid(): self.total_bounds = self.total_bounds.united(path.boundingRect()) pen = QPen(Qt.blue, 1) brush = QBrush(Qt.cyan) self.scene.addPath(path, pen, brush) data_source = None except Exception as e: print(f"加载失败: {str(e)}") self.current_vector_path = file_path # 新增这一行 data_source = None def geometry_to_qpainterpath(self, geom): path = QPainterPath() if geom.GetGeometryType() == ogr.wkbPolygon: for ring in range(geom.GetGeometryCount()): linear_ring = geom.GetGeometryRef(ring) points = linear_ring.GetPoints() if points: path.moveTo(points[0][0], points[0][1]) for p in points[1:]: path.lineTo(p[0], p[1]) path.closeSubpath() elif geom.GetGeometryType() == ogr.wkbLineString: points = geom.GetPoints() if points: path.moveTo(points[0][0], points[0][1]) for p in points[1:]: path.lineTo(p[0], p[1]) elif geom.GetGeometryType() == ogr.wkbPoint: x, y = geom.GetX(), geom.GetY() path.addEllipse(x - 2, y - 2, 4, 4) return path def auto_zoom(self): """自动缩放视图到数据范围并放大2倍""" if not self.total_bounds.isValid(): return # 设置场景边界 self.scene.setSceneRect(self.total_bounds) # 获取视图可视区域 view_rect = self.graphicsView.viewport().rect() # 计算缩放比例(自动适应 + 2倍放大) transform = QTransform() transform.scale(2, 2) # 先放大2倍 # 应用缩放并居中 self.graphicsView.setTransform(transform) self.graphicsView.fitInView(self.total_bounds, Qt.KeepAspectRatio) # 新增缓冲区分析方法 def buffer_analysis(self): """执行缓冲区分析""" if not hasattr(self, 'current_vector_path'): QMessageBox.warning(self, "警告", "请先打开矢量数据文件!") return # 获取缓冲距离 distance, ok = QInputDialog.getDouble( self, "缓冲区分析", "输入缓冲距离(单位与数据坐标系一致):", 0.0, 0 ) if not ok: return try: # 重新打开数据源获取几何 data_source = ogr.Open(self.current_vector_path, 0) layer = data_source.GetLayer(0) # 创建缓冲区路径 buffer_path = QPainterPath() pen = QPen(Qt.red, 2, Qt.DashLine) brush = QBrush(QColor(255, 0, 0, 50)) # 半透明红色填充 for feature in layer: geom = feature.GetGeometryRef() buffer_geom = geom.Buffer(distance) path = self.geometry_to_qpainterpath(buffer_geom) buffer_path.addPath(path) # 添加到场景 self.scene.addPath(buffer_path, pen, brush) # 更新视图边界 if buffer_path.boundingRect().isValid(): self.total_bounds = self.total_bounds.united(buffer_path.boundingRect()) self.auto_zoom() data_source = None except Exception as e: QMessageBox.critical(self, "错误", f"缓冲区分析失败: {str(e)}") def load_vector_data(self, file_path): self.scene.clear() self.total_bounds = QRectF() try: data_source = ogr.Open(file_path, 0) layer = data_source.GetLayer(0) # 获取字段定义 layer_defn = layer.GetLayerDefn() field_names = [layer_defn.GetFieldDefn(i).GetName() for i in range(layer_defn.GetFieldCount())] for feature in layer: geom = feature.GetGeometryRef() path = self.geometry_to_qpainterpath(geom) # 创建属性字典 attributes = { "FID": feature.GetFID(), **{name: feature.GetField(name) for name in field_names} } # 使用自定义图形项 item = FeatureItem(path, attributes) item.setPen(QPen(Qt.blue, 1)) item.setBrush(QBrush(Qt.cyan)) self.scene.addItem(item) if path.boundingRect().isValid(): self.total_bounds = self.total_bounds.united(path.boundingRect()) data_source = None except Exception as e: print(f"加载失败: {str(e)}") self.current_vector_path = file_path data_source = None # 新增属性查询方法 def enable_query_mode(self): """启用属性查询模式""" self.is_query_mode = not self.is_query_mode self.actionQuery_Attribute.setText("退出查询" if self.is_query_mode else "属性查询") self.graphicsView.setCursor(Qt.CrossCursor if self.is_query_mode else Qt.ArrowCursor) # 新增鼠标事件处理 def mousePressEvent(self, event): if self.is_query_mode and event.button() == Qt.LeftButton: scene_pos = self.graphicsView.mapToScene(event.pos()) items = self.scene.items(scene_pos, Qt.IntersectsItemShape, Qt.DescendingOrder) for item in items: if isinstance(item, FeatureItem): # 构建属性信息字符串 info = "\n".join([f"{k}: {v}" for k, v in item.attributes.items()]) QMessageBox.information(self, "要素属性", info) return super().mousePressEvent(event) def draw_centroids(self): """独立质心绘制功能""" if not hasattr(self, 'current_vector_path'): QMessageBox.warning(self, "警告", "请先打开矢量数据文件!") return # 清除已有质心 for item in self.centroid_items: self.scene.removeItem(item) self.centroid_items.clear() try: data_source = ogr.Open(self.current_vector_path, 0) layer = data_source.GetLayer(0) for feature in layer: geom = feature.GetGeometryRef() centroid = geom.Centroid() if centroid: # 创建质心图形项 path = QPainterPath() path.addEllipse( QRectF( centroid.GetX() - 0.3, # 修改为0.3像素半径 centroid.GetY() - 0.3, 0.6, 0.6 # 直径0.6像素 ) ) item = self.scene.addPath( path, QPen(Qt.red, 0.1), QBrush(Qt.red) ) self.centroid_items.append(item) data_source = None self.auto_zoom() except Exception as e: QMessageBox.critical(self, "错误", f"质心绘制失败: {str(e)}") # 新增空间查询模式切换方法 def enable_spatial_query_mode(self): """启用空间查询模式""" self.is_spatial_query_mode = not self.is_spatial_query_mode self.actionSpatial_Query.setText("退出空间查询" if self.is_spatial_query_mode else "空间查询") self.graphicsView.setCursor(Qt.CrossCursor if self.is_spatial_query_mode else Qt.ArrowCursor) if not self.is_spatial_query_mode: self.clear_spatial_query_results() # 新增空间查询处理方法 def mousePressEvent(self, event): if self.is_spatial_query_mode and event.button() == Qt.LeftButton: scene_pos = self.graphicsView.mapToScene(event.pos()) items = self.scene.items(scene_pos, Qt.IntersectsItemShape, Qt.DescendingOrder) for item in items: if isinstance(item, FeatureItem): # 获取空间关系选择 relations = ["相交", "包含", "被包含", "接触", "重叠"] relation, ok = QInputDialog.getItem( self, "空间关系选择", "请选择空间关系:", relations, 0, False ) if not ok: return # 执行空间查询 self.perform_spatial_query(item, relation) return super().mousePressEvent(event) # 新增空间查询核心方法 def perform_spatial_query(self, source_item, relation): """执行空间查询并高亮结果""" self.clear_spatial_query_results() try: # 获取源要素几何 source_geom = self.item_to_geometry(source_item) if not source_geom: return # 获取所有要素 all_items = [item for item in self.scene.items() if isinstance(item, FeatureItem)] # 遍历检查空间关系 for target_item in all_items: target_geom = self.item_to_geometry(target_item) if not target_geom: continue # 执行空间关系判断 if relation == "相交" and source_geom.Intersects(target_geom): self.highlight_item(target_item) elif relation == "包含" and source_geom.Contains(target_geom): self.highlight_item(target_item) elif relation == "被包含" and target_geom.Contains(source_geom): self.highlight_item(target_item) elif relation == "接触" and source_geom.Touches(target_geom): self.highlight_item(target_item) elif relation == "重叠" and source_geom.Overlaps(target_geom): self.highlight_item(target_item) except Exception as e: QMessageBox.critical(self, "错误", f"空间查询失败: {str(e)}") # 新增辅助方法 def item_to_geometry(self, item): """将图形项转换为OGR几何对象""" path = item.path() elements = path.toSubpathPolygons(QTransform()) if not elements: return None # 创建多边形几何 geom = ogr.Geometry(ogr.wkbPolygon) ring = ogr.Geometry(ogr.wkbLinearRing) for point in elements[0]: ring.AddPoint(point.x(), point.y()) ring.CloseRings() geom.AddGeometry(ring) return geom def highlight_item(self, item): """高亮显示查询结果""" original_pen = item.pen() highlight_pen = QPen(Qt.yellow, 3) item.setPen(highlight_pen) self.spatial_query_results.append((item, original_pen)) def clear_spatial_query_results(self): """清除查询结果高亮""" for item, original_pen in self.spatial_query_results: item.setPen(original_pen) self.spatial_query_results.clear() def open_raster_data(self): """打开栅格数据文件""" file_path, _ = QFileDialog.getOpenFileName( self, "打开栅格文件", "", "GeoTIFF (*.tif);;JPEG (*.jpg *.jpeg);;PNG (*.png);;All Files (*)" ) if file_path: try: self.load_raster_data(file_path) self.auto_zoom() except Exception as e: QMessageBox.critical(self, "错误", f"加载栅格失败: {str(e)}") def load_raster_data(self, file_path): """加载栅格数据到视图""" # 打开栅格文件(需要用户修改路径的部分) dataset = gdal.Open(file_path) # 相对路径示例:"./data/raster.tif" # 读取第一个波段 band = dataset.GetRasterBand(1) width = dataset.RasterXSize height = dataset.RasterYSize # 转换为numpy数组 data = band.ReadAsArray() # 创建QImage(注意数据类型转换) if data.dtype == np.uint8: format = QImage.Format.Format_Grayscale8 else: format = QImage.Format.Format_ARGB32 q_img = QImage(data.tobytes(), width, height, format) # 创建像素图项 pixmap = QPixmap.fromImage(q_img) raster_item = self.scene.addPixmap(pixmap) # 处理地理坐标(如果存在) geotransform = dataset.GetGeoTransform() if geotransform: # 计算四个角的坐标 x_origin = geotransform[0] y_origin = geotransform[3] pixel_width = geotransform[1] pixel_height = geotransform[5] # 更新场景边界 x_min = x_origin x_max = x_origin + pixel_width * width y_min = y_origin + pixel_height * height y_max = y_origin self.total_bounds = QRectF( QPointF(x_min, y_min), QPointF(x_max, y_max) ) dataset = None # 关闭数据集 def open_band_combination(self): if not hasattr(self, 'current_raster_path'): QMessageBox.warning(self, "警告", "请先打开栅格数据文件!") return # 复用open_raster_data的逻辑 self.open_raster_data() def open_raster_data(self): file_path, _ = QFileDialog.getOpenFileName( self, "打开栅格文件", "", "GeoTIFF (*.tif);;JPEG (*.jpg *.jpeg);;PNG (*.png);;All Files (*)" ) if file_path: try: dataset = gdal.Open(file_path) num_bands = dataset.RasterCount # 获取用户输入的波段组合 red_band, ok1 = QInputDialog.getInt( self, "波段选择", f"红通道波段号 (1-{num_bands}):", 1, 1, num_bands ) green_band, ok2 = QInputDialog.getInt( self, "波段选择", f"绿通道波段号 (1-{num_bands}):", min(2, num_bands), 1, num_bands ) blue_band, ok3 = QInputDialog.getInt( self, "波段选择", f"蓝通道波段号 (1-{num_bands}):", min(3, num_bands), 1, num_bands ) if not (ok1 and ok2 and ok3): return self.load_raster_data(file_path, red_band, green_band, blue_band) self.auto_zoom() self.current_raster_path = file_path # 新增存储当前路径 except Exception as e: QMessageBox.critical(self, "错误", f"加载栅格失败: {str(e)}") def load_raster_data(self, file_path, red_band=1, green_band=2, blue_band=3): """加载栅格数据到视图(支持波段组合)""" dataset = gdal.Open(file_path) width = dataset.RasterXSize height = dataset.RasterYSize # 读取三个波段数据 def read_band(band_num): band = dataset.GetRasterBand(band_num) data = band.ReadAsArray() # 自动拉伸到0-255范围 data_min = data.min() data_max = data.max() return np.clip(((data - data_min) / (data_max - data_min) * 255), 0, 255).astype(np.uint8) # 合并波段 rgb_array = np.dstack([ read_band(red_band), read_band(green_band), read_band(blue_band) ]) # 创建QImage q_img = QImage( rgb_array.data, width, height, 3 * width, # 每像素3字节(RGB) QImage.Format.Format_RGB888 ) # 创建像素图项 pixmap = QPixmap.fromImage(q_img) self.scene.addPixmap(pixmap) # 处理地理坐标(保持原有逻辑) geotransform = dataset.GetGeoTransform() if geotransform: x_origin = geotransform[0] y_origin = geotransform[3] pixel_width = geotransform[1] pixel_height = geotransform[5] x_min = x_origin x_max = x_origin + pixel_width * width y_min = y_origin + pixel_height * height # 计算下边界 y_max = y_origin # 上边界 # 确保坐标顺序正确 if x_min > x_max: x_min, x_max = x_max, x_min if y_min > y_max: y_min, y_max = y_max, y_min self.total_bounds = QRectF(QPointF(x_min, y_min), QPointF(x_max, y_max)) dataset = None # 新增栅格裁剪方法(必须缩进在类内部) def clip_raster(self): """执行栅格裁剪功能""" if not hasattr(self, 'current_raster_path'): QMessageBox.warning(self, "警告", "请先打开栅格数据文件!") return # 选择裁剪矢量文件 vector_path, _ = QFileDialog.getOpenFileName( self, "选择裁剪区域文件", "", "Shapefile (*.shp);;GeoJSON (*.geojson);;All Files (*)" ) if not vector_path: return try: # 获取原始栅格信息 src_ds = gdal.Open(self.current_raster_path) geotransform = src_ds.GetGeoTransform() proj = src_ds.GetProjection() # 获取矢量范围 vector_ds = ogr.Open(vector_path) layer = vector_ds.GetLayer() feature = layer.GetNextFeature() geom = feature.GetGeometryRef() x_min, x_max, y_min, y_max = geom.GetEnvelope() # 创建临时裁剪结果文件 import os # 确保导入os模块 output_path = os.path.splitext(self.current_raster_path)[0] + "_clipped.tif" # 执行裁剪操作 options = gdal.WarpOptions( format='GTiff', outputBounds=[x_min, y_min, x_max, y_max], dstSRS=proj ) gdal.Warp(output_path, src_ds, options=options) # 加载裁剪结果 self.load_raster_data(output_path) self.auto_zoom() # 清理资源 src_ds = None vector_ds = None except Exception as e: QMessageBox.critical(self, "错误", f"栅格裁剪失败: {str(e)}") # 新增波段运算方法 def band_calculation(self): """执行波段运算(示例为NDVI计算)""" if not hasattr(self, 'current_raster_path'): QMessageBox.warning(self, "警告", "请先打开栅格数据文件!") return try: # 获取用户输入参数 red_band, ok1 = QInputDialog.getInt( self, "波段选择", "输入红波段编号 (1-based):", 1, 1, 100 ) nir_band, ok2 = QInputDialog.getInt( self, "波段选择", "输入近红外波段编号 (1-based):", 4, 1, 100 ) if not (ok1 and ok2): return # 读取栅格数据 dataset = gdal.Open(self.current_raster_path) red = dataset.GetRasterBand(red_band).ReadAsArray() nir = dataset.GetRasterBand(nir_band).ReadAsArray() # 执行NDVI计算 ndvi = np.where( (nir + red) == 0, 0, (nir - red) / (nir + red) ).astype(np.float32) # 创建输出文件 output_path, _ = QFileDialog.getSaveFileName( self, "保存结果", "", "GeoTIFF (*.tif)" ) if not output_path: return # 写入结果 driver = gdal.GetDriverByName('GTiff') out_ds = driver.Create( output_path, dataset.RasterXSize, dataset.RasterYSize, 1, gdal.GDT_Float32 ) out_ds.SetGeoTransform(dataset.GetGeoTransform()) out_ds.SetProjection(dataset.GetProjection()) out_ds.GetRasterBand(1).WriteArray(ndvi) out_ds.FlushCache() # 清理资源 dataset = None out_ds = None QMessageBox.information(self, "成功", f"NDVI计算结果已保存至:\n{output_path}") except Exception as e: QMessageBox.critical(self, "错误", f"波段运算失败: {str(e)}") # 新增保存方法 def save_vector_data(self): """保存当前场景中的矢量数据""" if not hasattr(self, 'current_vector_path'): QMessageBox.warning(self, "警告", "没有可保存的矢量数据!") return # 获取保存路径 file_path, _ = QFileDialog.getSaveFileName( self, "保存矢量文件", "", "Shapefile (*.shp);;GeoJSON (*.geojson);;All Files (*)" ) if not file_path: return try: # 获取原始数据源信息 src_ds = ogr.Open(self.current_vector_path) src_layer = src_ds.GetLayer(0) src_defn = src_layer.GetLayerDefn() # 创建目标数据源 driver = ogr.GetDriverByName("ESRI Shapefile" if file_path.endswith(".shp") else "GeoJSON") if os.path.exists(file_path): driver.DeleteDataSource(file_path) dst_ds = driver.CreateDataSource(file_path) # 创建图层(保持与原始数据相同的坐标系) dst_layer = dst_ds.CreateLayer( "saved_features", srs=src_layer.GetSpatialRef(), geom_type=ogr.wkbPolygon ) # 复制字段定义 for i in range(src_defn.GetFieldCount()): field_defn = src_defn.GetFieldDefn(i) dst_layer.CreateField(field_defn) # 遍历场景中的所有要素项 for item in self.scene.items(): if isinstance(item, FeatureItem): # 创建新要素 feature = ogr.Feature(dst_layer.GetLayerDefn()) # 复制属性 for key, value in item.attributes.items(): if key == "FID": continue # FID通常自动生成 if feature.GetFieldIndex(key) != -1: feature.SetField(key, str(value)) # 转换几何 geom = self.item_to_geometry(item) if geom: feature.SetGeometry(geom) dst_layer.CreateFeature(feature) feature = None # 添加缓冲区要素(如果存在) self.save_additional_features(dst_layer, "buffer") # 添加质心要素(如果存在) self.save_additional_features(dst_layer, "centroid") dst_ds = None src_ds = None QMessageBox.information(self, "成功", f"数据已保存至:\n{file_path}") except Exception as e: QMessageBox.critical(self, "错误", f"保存失败: {str(e)}") # 新增辅助保存方法 def save_additional_features(self, layer, feature_type): """保存附加要素(缓冲区/质心)""" items = [] if feature_type == "buffer": items = [item for item in self.scene.items() if item.pen().style() == Qt.DashLine and item.pen().color() == Qt.red] elif feature_type == "centroid": items = self.centroid_items for item in items: geom = self.item_to_geometry(item) if geom: feature = ogr.Feature(layer.GetLayerDefn()) feature.SetGeometry(geom) # 添加类型标识字段 feature.SetField("FEATURE_TYPE", feature_type.upper()) layer.CreateFeature(feature) feature = None # 在item_to_geometry方法中添加点要素支持 def item_to_geometry(self, item): """增强版几何转换(支持点要素)""" path = item.path() elements = path.toSubpathPolygons(QTransform()) if not elements: # 处理点要素 if isinstance(item, QGraphicsPathItem): path = item.path() if path.elementCount() == 1 and path.elementAt(0).isMoveTo(): pt = path.elementAt(0) geom = ogr.Geometry(ogr.wkbPoint) geom.AddPoint(pt.x, pt.y) return geom return None # 原有多边形处理逻辑 geom = ogr.Geometry(ogr.wkbPolygon) ring = ogr.Geometry(ogr.wkbLinearRing) for point in elements[0]: ring.AddPoint(point.x(), point.y()) ring.CloseRings() geom.AddGeometry(ring) return geom # 新增栅格保存方法 def save_raster_data(self): """保存当前显示的栅格数据""" if not hasattr(self, 'current_raster_path'): QMessageBox.warning(self, "警告", "没有可保存的栅格数据!") return # 获取保存路径和格式 file_path, selected_filter = QFileDialog.getSaveFileName( self, "保存栅格文件", "", "GeoTIFF (*.tif);;JPEG (*.jpg *.jpeg);;PNG (*.png);;All Files (*)" ) if not file_path: return try: # 获取当前显示的栅格数据 dataset = gdal.Open(self.current_raster_path) band_count = dataset.RasterCount width = dataset.RasterXSize height = dataset.RasterYSize # 创建输出数据集 driver = gdal.GetDriverByName(self.get_driver_name(file_path)) out_ds = driver.Create( file_path, width, height, band_count, self.get_gdal_datatype(dataset) ) # 复制地理参考信息 out_ds.SetGeoTransform(dataset.GetGeoTransform()) out_ds.SetProjection(dataset.GetProjection()) # 复制所有波段数据 for band_num in range(1, band_count + 1): in_band = dataset.GetRasterBand(band_num) out_band = out_ds.GetRasterBand(band_num) out_band.WriteArray(in_band.ReadAsArray()) out_band.FlushCache() # 清理资源 out_ds = None dataset = None QMessageBox.information(self, "成功", f"栅格数据已保存至:\n{file_path}") except Exception as e: QMessageBox.critical(self, "错误", f"栅格保存失败: {str(e)}") # 新增辅助方法 def get_driver_name(self, file_path): """根据文件扩展名获取GDAL驱动名称""" ext = os.path.splitext(file_path)[1].lower() return { '.tif': 'GTiff', '.jpg': 'JPEG', '.jpeg': 'JPEG', '.png': 'PNG' }.get(ext, 'GTiff') def get_gdal_datatype(self, dataset): """获取与原始数据集匹配的GDAL数据类型""" sample_band = dataset.GetRasterBand(1) return { gdal.GDT_Byte: gdal.GDT_Byte, gdal.GDT_UInt16: gdal.GDT_UInt16, gdal.GDT_Int16: gdal.GDT_Int16, gdal.GDT_UInt32: gdal.GDT_UInt32, gdal.GDT_Int32: gdal.GDT_Int32, gdal.GDT_Float32: gdal.GDT_Float32, gdal.GDT_Float64: gdal.GDT_Float64 }.get(sample_band.DataType, gdal.GDT_Float32)
最新发布
07-18
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值