GeoTools入门（七）-- 操作栅格数据

亦凉 2022-10-10 15:45 1162阅读 0赞

## 1. 概述 ##

在前面的章节，我们主要讲述了GeoTools对于矢量(主要是Shape)数据的操作。在地理信息系统的世界里，还有一类很重要的数据类型，那便是栅格数据。虽然GeoTools对于栅格数据的支持并没有gdal强大，但既然他作为GeoTools的一部分，我们还是又必要了解一下它。

## 2. GridCoverage ##

GeoTools对于栅格数据的支持主要又由GridCoverage实现的。作为一名程序员，我们习惯于处理诸如JPEG、GIF或者PNG等格式的栅格数据。在地理空间方面，有一个Coverage个概念，他是空间定位要素的集合。非正式地，我们将地图和Coverage视为等同，当然，这是相对于地理意思而非编程思维上。

GridCoverage是Coverage的一个特例，他要求栅格要以矩形的形式填充到Coverage的区域中。在我们的java代码中，我们可以使用位图图形作为GridCoverage的后台数据结构，并使用其他元素记录特定坐标系中的空间边界。

这里有许多种GridCoverage的文件格式，最常用的如一下三类。

*  World+ 图像
    
    一种普通的图像格式，如jpeg或png，它有一个sidecar 文件来描述它的位置，还有一个`prjsidecar`文件来定义地图投影，就像shapefile使用的那样。请注意，尽管jpeg格式由于下载量小而很常见；运行时的性能非常糟糕，因为整个图像都需要读入内存。TIFF等格式没有此限制。
 *  GeoTiff
    
    具有存储在图像元数据字段中的地理空间信息的普通tiff图像。这通常在安全的前提下表现出良好的性能；特别是如果它已经准备了一个内部覆盖（可用于缩小）或内部平铺（允许快速平移时放大）。当计算机的磁盘速度比CPU快时，性能最佳。
 *  JPEG2000
    
    他jpeg的扩展，使用小波压缩来处理大量图像。文件格式还支持可用于存储地理空间信息的元数据字段。当您CPU速度比磁盘访问更快的情况下，这种格式的性能最好。

针对上面三种数据个数的描述，我们知道，第一种数据格式基本是被摒弃的那种。GeoTiff和JPEG2000的区别是要看IO性能VS CPU性能。

## 3. GridCoverage对栅格数据的操作 ##

在本篇文章中，我们主要关注GridCoverage对于栅格数据的加载，波段提取，显示样式的创建。

### 3.1 栅格数据的加载 ###

AbstractGridFormat format = GridFormatFinder.findFormat(rasterFile);
    // this is a bit hacky but does make more geotiffs work
    Hints hints = new Hints();
    if (format instanceof GeoTiffFormat) { 
        hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
    }
    reader = format.getReader(rasterFile, hints);

### 3.2 波段提取 ###

GridCoverage2D cov = reader.read(null);
    int numBands = cov.getNumSampleDimensions();
    for (int i = 0; i < numBands; i++) { 
        GridSampleDimension dim = cov.getSampleDimension(i);
        //这里依次输出 RED_BAND GREEN_BANK BLUE_BANK
        System.out.println(dim.getDescription().toString());
    }

### 3.3 样式创建 ###

#### 3.3.1 创建灰度的样式 ####

private Style createGreyscaleStyle(int band) { 
        StyleFactory sf = CommonFactoryFinder.getStyleFactory();
        FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
        ContrastEnhancement ce = sf.contrastEnhancement(ff.literal(1.0), ContrastMethod.NORMALIZE);
        SelectedChannelType sct = sf.createSelectedChannelType(String.valueOf(band), ce);
        RasterSymbolizer sym = sf.getDefaultRasterSymbolizer();
        ChannelSelection sel = sf.channelSelection(sct);
        sym.setChannelSelection(sel);
        return SLD.wrapSymbolizers(sym);
    }

#### 3.3.2 创建RGB样式 ####

StyleFactory sf = CommonFactoryFinder.getStyleFactory();
    FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
    SelectedChannelType[] sct = new SelectedChannelType[cov.getNumSampleDimensions()];
    ContrastEnhancement ce = sf.contrastEnhancement(ff.literal(1.0), ContrastMethod.NORMALIZE);
    for (int i = 0; i < 3; i++) { 
        sct[i] = sf.createSelectedChannelType(String.valueOf(channelNum[i]), ce);
    }
    RasterSymbolizer sym = sf.getDefaultRasterSymbolizer();
    ChannelSelection sel = sf.channelSelection(sct[RED], sct[GREEN], sct[BLUE]);
    sym.setChannelSelection(sel);
    return SLD.wrapSymbolizers(sym);

通过上面的例子，我们可以简单画个时序图及流程图方便大家记忆：

CommonFactoryFinderStyleFactoryFilterFactory2ContrastEnhancementSelectedChannelTypeRasterSymbolizerChannelSelectionSLDStyle1.getStyleFactory2.getFilterFactory23.contrastEnhancement(CommonFactoryFinder)4.createSelectedChannelType5.getDefaultRasterSymbolizer6.channelSelection(SelectedChannelType)7.setChannelSelection(ChannelSelection)7.wrapSymbolizers(RasterSymbolizer)CommonFactoryFinderStyleFactoryFilterFactory2ContrastEnhancementSelectedChannelTypeRasterSymbolizerChannelSelectionSLDStyle

## 4. 示例 ##

前面的示例中，我们查看了如何读取和显示Shape文件。对于ImageLab.java，我们将通过显示一个三波段的全球卫星图像，并将其与shapefile中的国家边界进行叠加，将光栅数据添加到混合中。

### 4.1 pom引入 ###

<dependencies>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-shapefile</artifactId>
            <version>${geotools.version}</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-swing</artifactId>
            <version>${geotools.version}</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-epsg-hsql</artifactId>
            <version>${geotools.version}</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-geotiff</artifactId>
            <version>${geotools.version}</version>
        </dependency>
        <dependency>
            <groupId>org.geotools</groupId>
            <artifactId>gt-image</artifactId>
            <version>${geotools.version}</version>
        </dependency>
    </dependencies>
    <repositories>
        <repository>
            <id>osgeo</id>
            <name>OSGeo Release Repository</name>
            <url>https://repo.osgeo.org/repository/release/</url>
            <snapshots><enabled>false</enabled></snapshots>
            <releases><enabled>true</enabled></releases>
        </repository>
        <repository>
            <id>osgeo-snapshot</id>
            <name>OSGeo Snapshot Repository</name>
            <url>https://repo.osgeo.org/repository/snapshot/</url>
            <snapshots><enabled>true</enabled></snapshots>
            <releases><enabled>false</enabled></releases>
        </repository>
    </repositories>

### 4.2 创建ImageLab类 ###

public class ImageLab { 
        private StyleFactory sf = CommonFactoryFinder.getStyleFactory();
        private FilterFactory2 ff = CommonFactoryFinder.getFilterFactory2();
        private JMapFrame frame;
        private GridCoverage2DReader reader;
    
        public static void main(String[] args) throws Exception { 
            ImageLab me = new ImageLab();
            me.getLayersAndDisplay();
        }
    }

### 4.3 选择Image和shape文件 ###

private void getLayersAndDisplay() throws Exception { 
         List<Parameter<?>> list = new ArrayList<>();
         list.add(
             new Parameter<>(
                 "image",
                 File.class,
                 "Image",
                 "GeoTiff or World+Image to display as basemap",
                 new KVP(Parameter.EXT, "tif", Parameter.EXT, "jpg")));
         list.add(
             new Parameter<>(
                 "shape",
                 File.class,
                 "Shapefile",
                 "Shapefile contents to display",
                 new KVP(Parameter.EXT, "shp")));
    
         JParameterListWizard wizard =
             new JParameterListWizard("Image Lab", "Fill in the following layers", list);
         int finish = wizard.showModalDialog();
    
         if (finish != JWizard.FINISH) { 
             System.exit(0);
         }
         File imageFile = (File) wizard.getConnectionParameters().get("image");
         File shapeFile = (File) wizard.getConnectionParameters().get("shape");
         displayLayers(imageFile, shapeFile);
     }

### 4.4 显示地图 ###

private void displayLayers(File rasterFile, File shpFile) throws Exception { 
        AbstractGridFormat format = GridFormatFinder.findFormat(rasterFile);
        Hints hints = new Hints();
        if (format instanceof GeoTiffFormat) { 
            hints = new Hints(Hints.FORCE_LONGITUDE_FIRST_AXIS_ORDER, Boolean.TRUE);
        }
        reader = format.getReader(rasterFile, hints);
    
        // 创建一个灰度样式
        Style rasterStyle = createGreyscaleStyle(1);
    
        // 连接shape图层
        FileDataStore dataStore = FileDataStoreFinder.getDataStore(shpFile);
        SimpleFeatureSource shapefileSource = dataStore.getFeatureSource();
    
        //创建一个黄色边界且没有填充的简单样式
        Style shpStyle = SLD.createPolygonStyle(Color.YELLOW, null, 0.0f);
    
        // 创建一个两个图层的地图
        final MapContent map = new MapContent();
        map.setTitle("ImageLab");
    
        Layer rasterLayer = new GridReaderLayer(reader, rasterStyle);
        map.addLayer(rasterLayer);
    
        Layer shpLayer = new FeatureLayer(shapefileSource, shpStyle);
        map.addLayer(shpLayer);
    
        //设置地图界面的样式
        frame = new JMapFrame(map);
        frame.setSize(800, 600);
        frame.enableStatusBar(true);
        frame.enableToolBar(true);
    
        JMenuBar menuBar = new JMenuBar();
        frame.setJMenuBar(menuBar);
        JMenu menu = new JMenu("Raster");
        menuBar.add(menu);
    
        menu.add(
            new SafeAction("Grayscale display") { 
                public void action(ActionEvent e) throws Throwable { 
                    Style style = createGreyscaleStyle();
                    if (style != null) { 
                        ((StyleLayer) map.layers().get(0)).setStyle(style);
                        frame.repaint();
                    }
                }
            });
    
        menu.add(
            new SafeAction("RGB display") { 
                public void action(ActionEvent e) throws Throwable { 
                    Style style = createRGBStyle();
                    if (style != null) { 
                        ((StyleLayer) map.layers().get(0)).setStyle(style);
                        frame.repaint();
                    }
                }
            });
        frame.setVisible(true);
    }

#### 4.4.1 创建选择灰度的对话框 ####

private Style createGreyscaleStyle() { 
        GridCoverage2D cov = null;
        try { 
            cov = reader.read(null);
        } catch (IOException giveUp) { 
            throw new RuntimeException(giveUp);
        }
        int numBands = cov.getNumSampleDimensions();
        Integer[] bandNumbers = new Integer[numBands];
        for (int i = 0; i < numBands; i++) { 
            bandNumbers[i] = i + 1;
        }
        Object selection =
            JOptionPane.showInputDialog(
            frame,
            "Band to use for greyscale display",
            "Select an image band",
            JOptionPane.QUESTION_MESSAGE,
            null,
            bandNumbers,
            1);
        if (selection != null) { 
            int band = ((Number) selection).intValue();
            return createGreyscaleStyle(band);
        }
        return null;
    }

#### 4.4.2 创建灰度样式 ####

private Style createGreyscaleStyle(int band) { 
        ContrastEnhancement ce = sf.contrastEnhancement(ff.literal(1.0), ContrastMethod.NORMALIZE);
        SelectedChannelType sct = sf.createSelectedChannelType(String.valueOf(band), ce);
    
        RasterSymbolizer sym = sf.getDefaultRasterSymbolizer();
        ChannelSelection sel = sf.channelSelection(sct);
        sym.setChannelSelection(sel);
    
        return SLD.wrapSymbolizers(sym);
    }

### 4.5 创建RGB样式 ###

private Style createRGBStyle() { 
        GridCoverage2D cov = null;
        try { 
            cov = reader.read(null);
        } catch (IOException giveUp) { 
            throw new RuntimeException(giveUp);
        }
        //我们需要至少三个波段的样式
        int numBands = cov.getNumSampleDimensions();
        if (numBands < 3) { 
            return null;
        }
        //获取各个波段的名称
        String[] sampleDimensionNames = new String[numBands];
        for (int i = 0; i < numBands; i++) { 
            GridSampleDimension dim = cov.getSampleDimension(i);
            sampleDimensionNames[i] = dim.getDescription().toString();
        }
        final int RED = 0, GREEN = 1, BLUE = 2;
        int[] channelNum = { -1, -1, -1};
        for (int i = 0; i < numBands; i++) { 
            String name = sampleDimensionNames[i].toLowerCase();
            if (name != null) { 
                if (name.matches("red.*")) { 
                    channelNum[RED] = i + 1;
                } else if (name.matches("green.*")) { 
                    channelNum[GREEN] = i + 1;
                } else if (name.matches("blue.*")) { 
                    channelNum[BLUE] = i + 1;
                }
            }
        }
        if (channelNum[RED] < 0 || channelNum[GREEN] < 0 || channelNum[BLUE] < 0) { 
            channelNum[RED] = 1;
            channelNum[GREEN] = 2;
            channelNum[BLUE] = 3;
        }
        //我们使用selected通道创建RasterSymbolizer样式
        SelectedChannelType[] sct = new SelectedChannelType[cov.getNumSampleDimensions()];
        ContrastEnhancement ce = sf.contrastEnhancement(ff.literal(1.0), ContrastMethod.NORMALIZE);
        for (int i = 0; i < 3; i++) { 
            sct[i] = sf.createSelectedChannelType(String.valueOf(channelNum[i]), ce);
        }
        RasterSymbolizer sym = sf.getDefaultRasterSymbolizer();
        ChannelSelection sel = sf.channelSelection(sct[RED], sct[GREEN], sct[BLUE]);
        sym.setChannelSelection(sel);
    
        return SLD.wrapSymbolizers(sym);
    }