时代互联资讯中心
为什么极值如此关键
做地理数据处理的同行多半都有过这样的经历：手里攥着一个动辄几百M的GeoTIFF文件，第一句话往往想问的是——这组数据到底覆盖了多大数值区间？GeoTIFF作为地理空间数据的主流载体，里面塞着人口、地形、夜间灯光这类数据，从遥感像片到地貌模型无所不包。地图绘制、资源勘查、环境监测等领域都离不开它。拿到手的第一步，很多人并不是急着建模，而是要摸清数据的分布底细，而最大值和最小值正是这条底线上最醒目的两个端点。
这个范围信息看着朴素，实际用处却不少。拿遥感来说，不同地物反射率、发射率高低不等，知道峰值谷值，就能对地物构成和数据本身的质量有个初步判断。地形方面就更直接了，高程的最大最小值相减，差不多就等于整个研究区的海拔落差，后面要算坡度、坡向，都绕不开它打底子。
要拿到这两个值，路子其实有两条：一条是图形界面，靠QGIS这类软件点开就能看；另一条是写代码，用开源库GeoTools自己算。两者结论应该对得上。下面先把界面版讲清楚，再重点说说怎么在程序里把这件事自动化。
第一招：用QGIS界面直接读数
如果你只是偶尔查一次，完全不必碰代码。QGIS加载tif之后，有好几个地方都能看到极值，顺便说一句，ArcMap、ArcGIS那套界面里也有类似入口。这里挑QGIS最常用的三种展示位置讲。
其一：属性信息标签页
把tif拖进QGIS，右键点开图层属性，切到属性信息标签页。界面右侧会把该波段的一些统计量直接摊给你看，最大、最小这些字段通常就在其中，翻一翻就能查到。
这里列出的是该波段自带的统计结果，打开就能看到数值。
STATISTICS_MAXIMUM=2095
STATISTICS_MEAN=355.64
STATISTICS_MINIMUM=-35
STATISTICS_STDDEV=297.91
STATISTICS_VALID_PERCENT=64.56
上面的记录就是波段统计区的输出样貌：最大值两千零九十五，最小值负三十五，中间还夹着均值、标准差和有效像素占比。字段含义一目了然。
其二：符号化标签页
平时给图层上颜色的时候会进到符号化标签页。这里有一段特殊值输入区，下拉框右侧会把当前波段的范围补齐给你看，最大、最小连着一起列着，改拉伸区间时顺手就用了。
顺着这个位置读到的最大值同样也是两千零九十五，最小值负三十五，和刚才属性信息里对得上号，说明数据是自洽的。
其三：直方图统计页
栅格图层比矢量多了个兄弟标签页——直方图。点进去会弹出一张柱状分布图，图下方同样标注着当前tif对应的最大、最小值。
这张图把数据分布摊得更直观，底部标出来的端点数值，依旧和前面两种方式一致。
第二招：用GeoTools把过程自动化
界面查值适合人看，但一旦要做批量处理、要嵌进pipeline，就得分交出代码了。GeoTools这个开源库大家都不陌生，扩展性强，跟各类GIS格式都能打交道，算法模块也够健壮。用它来算GeoTIFF的极值，既省了自己造轮子，效率和精度也更有保障。下面分享两套常用写法：一套是老老实实遍历像素，另一套是走统计分类的现成功能。
方式一：直接遍历像素
这套思路最直白，一层层读栅格坐标点，挨个过一遍。先把像元读进数组，再过滤掉无效值，把有用的点攒起来排序，头尾两个自然就是我们要的极值。思路简单到有点天真，小范围数据跑得挺欢，但碰到几百平方公里的大图，整个内存很容易就被吃空了。
下面是这套写法的骨架：
File file = new File("F:/vector_data/worldPop人口数据/worldpop/chn_ppp_2020_1km_Aggregated.tif");
GeoTiffReader reader = new GeoTiffReader(file, new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE));
GridCoverage2D coverage = reader.read(null);
// 重采样到目标坐标系
CoordinateReferenceSystem targetCRS = CRS.decode("EPSG:4326");
coverage = (GridCoverage2D) Operations.DEFAULT.resample(coverage, targetCRS);
RenderedImage image = coverage.getRenderedImage();
GridGeometry2D geometry = coverage.getGridGeometry();
// 整块读进内存
Raster raster = image.getData();
width = raster.getWidth();
int height = raster.getHeight();
double[] pixels = raster.getSamples(raster.getMinX(), raster.getMinY(), width, height, 0, (double[]) null);
// 只挑非零的坐标攒起来
List nonZeroCoords = new ArrayList<>();
for (int j = 0; j < height; j++) { for (int i = 0; i < width; i++) { if (pixels[j * width + i] > 0) {
nonZeroCoords.add(new GridCoordinates2D(i, j));
}
}
}
MathTransform2D gridToCRS = geometry.getGridToCRS2D();
if (!(gridToCRS instanceof AffineTransform2D)) {
throw new RuntimeException("仅支持仿射变换");
}
AffineTransform2D affine = (AffineTransform2D) gridToCRS;
BigDecimal total = new BigDecimal("0.0");
double t = 0;
// 把像素中心坐标反算成经纬度
for (GridCoordinates2D coord : nonZeroCoords) {
int i = coord.x;
int j = coord.y;
double[] src = new double[]{i + 0.5, j + 0.5};
double[] dest = new double[2];
affine.transform(src, 0, dest, 0, 1);
double longitude = dest[0];
double latitude = dest[1];
total = total.add(new BigDecimal(pixels[j * width + i]));
t += pixels[j * width + i];
}
这套写法的用武之地，主要是你想顺带拿到每个像元的坐标归属，比如人口密度、高程这类数据既要数值也要位置。但如果只在乎那个最大最小值，就显得杀鸡用牛刀了。这里给个实在的建议：大批量场景别硬刚，见文末那个统计类，能省掉大量体力活。顺带一提，要是你的算力环境有点紧，云主机怎么选、怎么配内存，去时代互联看看方案也是个省心的路子。
方式二：走CoverageClassStats统计
还记得QGIS属性信息里那段自带统计吗？GeoTools里就有对应的现成功能，不必自己遍历。借助CoverageClassStats，把最大、最小两个统计量扔进去，一把算完。代码長相这样：
package com.yelang.project.geotoolstiff;
import java.io.File;
import java.io.IOException;
import java.util.HashSet;
import java.util.Set;
import org.geotools.coverage.grid.GridCoverage2D;
import org.geotools.coverage.grid.io.AbstractGridCoverage2DReader;
import org.geotools.coverage.grid.io.AbstractGridFormat;
import org.geotools.coverage.grid.io.GridFormatFinder;
import org.geotools.gce.geotiff.GeoTiffFormat;
import org.geotools.gce.geotiff.GeoTiffReader;
import org.geotools.process.ProcessException;
import org.geotools.process.classify.ClassificationMethod;
import org.geotools.process.raster.CoverageClassStats;
import org.geotools.util.factory.Hints;
import org.jaitools.numeric.Statistic;
public class DemMinMaxStat {
public static void main(String[] args) throws ProcessException, IOException {
// 湖南地形 min：-35 max：2095
File file = new File("C:/BaiduDownload/湖南省_DEM_30m分辨率_NASA数据.tif");
// 全国人口 min：-99999.0 max：286414.71875
//File file = new File("F:/vector_data/worldPop人口数据/worldpop/chn_ppp_2020_1km_Aggregated.tif");
AbstractGridFormat format = GridFormatFinder.findFormat(file);
if (!(format instanceof GeoTiffFormat)) {
throw new RuntimeException(String.format("传入文件不是geoTif格式文件"));
}
Hints hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
AbstractGridCoverage2DReader reader = format.getReader(file, hints);
GridCoverage2D coverage = reader.read(null);
CoverageClassStats rasterProcess = new CoverageClassStats();
Set set = new HashSet(2);
set.add(Statistic.MAX); // 最大值
set.add(Statistic.MIN); // 最小值
CoverageClassStats.Results results = rasterProcess.execute(coverage, set, 0, 1,
ClassificationMethod.QUANTILE, null, null);
System.out.println("geometry = " + results);
double maxVal = results.value(0, Statistic.MAX);
double minVal = results.value(0, Statistic.MIN);
System.out.println("Minimum value: " + minVal);
System.out.println("Maximum value: " + maxVal);
}
}
在IDE里直接跑，控制台就会把结果吐出来。拿湖南那笔DEM数据来说，读出来的最大值正好两千零九十五，最小值负三十五；换成全国人口那张图，则是最小负九万九千九百九十九、最大接近二十八万六千四百。跟QGIS界面里看到的分毫不差，结果的可靠性这就立住了。
这里有个小细节值得留意：CLASSIFICATION Method这里用的是QUANTILE分位数拉伸，如果你想要的是MIN_MAX线性拉伸，把参数换一下就行。统计口径不同，结果可能会有出入，别到时候对不上号又找不到原因。
写在最后
回到开头那个问题——怎么拿到GeoTIFF的最大最小值？界面党用QGIS，点开属性、符号化或直方图任何一个标签页都能读数；代码党则可以在遍历像素和调用CoverageClassStats之间二选一。前者能顺带拿到空间坐标，适合小范围精细活；后者纯算统计量，跑大图更从容。两种路子最终要回到同一个数上，互相对账，才是靠谱的做法。
技术这东西，纸上得来终觉浅。建议大家拿手头的真实tif先跑一遍，尤其是数据量大的场景，更要亲自测一下内存占用。文章到这里就算讲完了，若有疏漏，欢迎留言拍砖，共同进步。
