GDAL+Java实现获取对应栅格影像经纬度对应的像素值 |
您所在的位置：网站首页 › 纬度和经度的表达方法 › GDAL+Java实现获取对应栅格影像经纬度对应的像素值 |
GDAL+Java实现获取对应栅格影像经纬度对应的像素值
|
从前面的GDAL系列博文中，可以指导GDAL可以将栅格影像文件读出为对应的多维数组，可以读出每一个像素格对应的像素值。但如何根据经纬度直接读取像素值呢？博主从查阅了网上的相关文档，发现有个人写的计算公式是错误的，用代码跑出来的结果都是错误的。于是自己查阅了相关文档，自己实现了一遍，跟大家分享一下。 图像坐标空间到地理参考坐标空间的转换可以用GDAL读取出每个影像文件的空间地理变换，地理变换是从图像坐标空间（行、列），也称为（像素、线）到地理参考坐标空间（投影或地理坐标）的仿射变换。可以用GetGeoTransform读取出来，由以下六个参数组成 GT(0) 左上像素左上角的x坐标。GT(1) w-e像素分辨率/像素宽度。GT(2) 行旋转（通常为零）。GT(3) 左上像素左上角的y坐标。GT(4) 列旋转（通常为零）。GT(5) n-s像素分辨率/像素高度（北上图像为负值）。从图像坐标空间到地理参考坐标空间的转换关系如下： X_geo = GT(0) + X_pixel * GT(1) + Y_line * GT(2) Y_geo = GT(3) + X_pixel * GT(4) + Y_line * GT(5)X_geo为实际X坐标，Y_geo为实际Y坐标，X_pixel代表栅格数据的第n列，Y_line代表栅格数据的第n行。 注意，上面的像素/线坐标是从左上角的（0.0,0.0）到右下角的（宽度单位为像素，高度单位为像素）。因此，左上角像素中心的像素/线位置为（0.5,0.5）。 如果是北向图像： GT(2) ， GT(4) 系数为零。 GT(1) ， GT(5) 是像素大小。 GT(0) ， GT(3) 位置是栅格左上角像素的左上角。 计算公式从前面的关系中，可以倒退出对应的计算关系，和计算二元一次方程是一样的。 X_pixel =(GT(5)*X_geo-GT(2)*Y_geo-GT(5)*GT(0)+GT(2)*GT(3))/(GT(5)*GT(1)-GT(2)*GT(4)) Y_line =(Y_geo-GT(3)-xPixel*GT(4))/GT(5)注意： GT(2) ， GT(4)不能作为被除数，这可能导致计算错误 代码实现 package com.meteorology.clock.common; import org.gdal.gdal.Band; import org.gdal.gdal.Dataset; import org.gdal.gdal.Driver; import org.gdal.gdal.gdal; import org.gdal.gdalconst.gdalconstConstants; import java.util.Arrays; public class RasterTool { public static void main(String[] args) { //指定经纬度 double Latitude = 28.850445; double longitude = 119.100341; String filepath = "E:\\Z_NAFP_C_BABJ_20230712180520_P_HRCLDAS_RT_BEZZ_0P01_HOR-WIND-2023071218.img"; //注册 gdal.AllRegister(); //打开文件获取数据集 Dataset dataset = gdal.Open(filepath, gdalconstConstants.GA_ReadOnly); if (dataset == null) { System.out.println("打开"+filepath+"失败"+gdal.GetLastErrorMsg()); System.exit(1); } //获取驱动 Driver driver = dataset.GetDriver(); //获取驱动信息 System.out.println("driver long name: " + driver.getLongName()); //获取栅格数量 int bandCount = dataset.getRasterCount(); System.out.println("RasterCount: " + bandCount); //构造仿射变换参数数组，并获取数据 double[] gt = new double[6]; dataset.GetGeoTransform(gt); System.out.println("仿射变换参数"+ Arrays.toString(gt)); //经纬度转换为栅格像素坐标 int[] ColRow=Coordinates2ColRow(gt,longitude,Latitude); //判断是否坐标超出范围 if(ColRow[0]=dataset.getRasterYSize()){ System.out.println(Arrays.toString(ColRow)+"坐标值超出栅格范围！"); return; } //遍历波段，获取该点对应的每个波段的值并打印到屏幕 for (int i = 0;i // GT(0) 左上像素左上角的x坐标。 // GT(1) w-e像素分辨率/像素宽度。 // GT(2) 行旋转（通常为零）。 // GT(3) 左上像素左上角的y坐标。 // GT(4) 列旋转（通常为零）。 // GT(5) n-s像素分辨率/像素高度（北上图像为负值） // 如果是北向图像： // GT(2) ， GT(4) 系数为零。 // GT(1) ， GT(5) 是像素大小。 // GT(0) ， GT(3) 位置是栅格左上角像素的左上角。 // X_geo = GT(0) + X_pixel * GT(1) + Y_line * GT(2) // Y_geo = GT(3) + X_pixel * GT(4) + Y_line * GT(5) int[] ints = new int[2]; //向下取整,如果向上取整会导致计算结果偏大，从而在后面读取到邻近像元的数据 double xPixel=(gt[5]*X-gt[2]*Y-gt[5]*gt[0]+gt[2]*gt[3])/(gt[5]*gt[1]-gt[2]*gt[4]); double yLine=(Y-gt[3]-xPixel*gt[4])/gt[5]; ints[0]= (int) Math.floor(xPixel); ints[1]=(int) Math.floor(yLine); System.out.println(ints[0]); System.out.println(ints[1]); return ints; } } 结果验证使用ArcGIS工具，打开栅格影像，并查询对应的像素值 我这里使用的测试栅格影像是WGS84坐标系的，所以这样的计算没有问题，要是其他坐标系的，需要进行坐标转换。 |
程序结果： 
从上面可以看出，结果一致，而且小数位比ArcGIS更多。【本文地址】
公司简介
联系我们
今日新闻 |
点击排行 |
|
推荐新闻 |
图片新闻 |
|
专题文章 |