基于GDAL和wgrib2处理grib2气象文件

1.流域边界shp文件处理

把每一个流域的wata.shp文件合并一个整体流域的shp文件，放到要裁剪的shp文件夹

2.GRIB2—>NC文件

	/**
     * grib转化为nc
     *
     * @param inputFile    输入的grib2文件
     * @param outTifFolder 输出文件夹
     * @return
     */
    public static String grib2ToNc(String inputFile, String outTifFolder) {
        File tempFile = new File(inputFile.trim());
        //Z_NWGD_C_BEKM_20220623061226_P_RFFC_SPCC-ER03_202206230800_02403.GRB2

        if (!tempFile.isFile()) {
            throw new CustomException(ExceptionMessageEnum.GRIB2_NOT_FOUND);
        }
        String inputFileName = tempFile.getName();
        //剔除出来nc文件的名称 2022062308.nc
        String ncName = inputFileName.split("_")[8].substring(0, 10);
        String outPutNcFileName = ncName + ".nc";

        //设置nc文件保存的全路径 C:\Users\Desktop\GRIB_OUT\2022062308\2022062308.nc
        StringBuilder outPutNcFile = new StringBuilder(outTifFolder).append(GdalConstants.SPLIT).append(outPutNcFileName);


        //创建tif文件的输出文件夹, 如果当前路径不存在, 就创建
        File file = new File(outTifFolder + GdalConstants.SPLIT + ncName);
        if (!file.exists()) {
            file.mkdir();
        }

        // wgrib2 C:\Users\Desktop\GRIB\Z_NWGD_C_BEKM_20220918060027_P_RFFC_SPCC-ER03_202209180800_02403.GRB2 -netcdf C:\Users\Desktop\GRIB_OUT\2022091808.nc
        String cmdGrb2ToNc = String.format("wgrib2 %s -netcdf %s", inputFile, outPutNcFile);
        log.info("执行的转出nc脚本: " + cmdGrb2ToNc);
        //执行命令行命令
        boolean exeFlag = ProcessUtil.execCommand(cmdGrb2ToNc);
        if (exeFlag) {
            log.info("{} 文件转出成功", outPutNcFile);
            return outPutNcFile.toString();
        } else {
            log.info("{} 文件转出失败", outPutNcFile);
            return null;
        }
    }

3.NC文件裁剪

/**
     * tif重采样方法
     *
     * @param inputNcFile  输入的ncFile
     * @param outTifFolder 输出tif的文件夹
     * @param shpFolder    用于裁剪的shp文件夹
     */
    private Boolean ncToReSamTif(String inputNcFile, String outTifFolder, String shpFolder) {
        try {
            //使用netCDF读取nc文件
            NetcdfFile ncDataset = NetcdfFile.open(inputNcFile);
            Variable apcpVar = ncDataset.findVariable(GdalConstants.APCP_SURFACE);

            double[] geoTransform = getGeoTransform(ncDataset, 1);

            int nLat = apcpVar.getShape()[1];
            int nLon = apcpVar.getShape()[2];

            gdal.AllRegister();
            Driver driver = gdal.GetDriverByName(GdalConstants.DRIVER_NAME);
            //创建tif文件,这里比较麻烦:
            //1.创建8-11,11-14,14-17的8个tif文件
            //2.创建8-14,14-20,20-02的4个tif文件
            //3.创建24小时的tif文件
            String[] nameSuffix = {"003", "006", "024"};

            for (String suffix : nameSuffix) {
                //此时代表一个时间片导出为一个tif文件
                if ("003".equals(suffix)) {
                    for (int i = 1; i <= 8; i++) {
                        float[][] apcpArray = readNcFile(apcpVar, i - 1, 1);
                        //创建tif文件名称 2022082808.003.1.tif 2022082808.003.2.tif
                        String tifName = makeFileName(outTifFolder, suffix, i);
                        genTif(outTifFolder, shpFolder, ncDataset, geoTransform, nLat, nLon, driver, apcpArray, tifName);
                    }
                } else if ("006".equals(suffix)) {
                    //此时代表两个时间片导出为一个tif文件 12 34 56 78
                    for (int i = 1; i <= 4; i++) {
                        float[][] apcpArray = readNcFile(apcpVar, (i - 1) * 2, 2);
                        //创建tif文件名称 2022082808.006.1.tif  2022082808.006.2.tif
                        String tifName = makeFileName(outTifFolder, suffix, i);
                        genTif(outTifFolder, shpFolder, ncDataset, geoTransform, nLat, nLon, driver, apcpArray, tifName);
                    }
                } else if ("024".equals(suffix)) {
                    //此时代表两个时间片导出为一个tif文件 12 34 56 78
                    for (int i = 1; i <= 1; i++) {
                        float[][] apcpArray = readNcFile(apcpVar, 0, 8);
                        //创建tif文件名称 2022082808.024.1.tif
                        String tifName = makeFileName(outTifFolder, suffix, i);
                        genTif(outTifFolder, shpFolder, ncDataset, geoTransform, nLat, nLon, driver, apcpArray, tifName);
                    }
                }
            }
            return Boolean.TRUE;
        } catch (Exception e) {
            e.printStackTrace();
        }
        return Boolean.FALSE;
    }

读取nc文件的方法：

 /**
     * 从气象nc数据中裁剪指定区域的数据
     *
     * @param apcpVar 雨量变量
     * @param start   从第几片开始读取
     * @param size    读取多少片
     * @return 读取到的结果
     */
    public static float[][] readNcFile(Variable apcpVar, int start, int size) {
        try {
            float[][] values = null;
            //本读取的是一个3维数据，时间、经度、维度，一下子把3维数据全部读出来会很大，时间维度是8层，所以通过遍历时间维度获取数据，i为时间维度的层数
            for (int i = 0; i < size; i++) {
                //origin 设置维度准备从哪个位置开始读取数据
                int[] origin = {start + i, 0, 0};
                //shape 和origin对应，设置读取多少点后结束
                int[] shape = apcpVar.getShape();
                shape[0] = 1;
                //去掉时间维度，变为二维
                float[][] temp = (float[][]) apcpVar.read(origin, shape).reduce(0).copyToNDJavaArray();
                if (values != null) {
                    for (int j = 0; j < shape[1]; j++) {
                        for (int k = 0; k < shape[2]; k++) {
                            values[j][k] += temp[j][k];
                        }
                    }
                } else {
                    values = temp;
                }

            }
            return values;
        } catch (IOException | InvalidRangeException e) {
            e.printStackTrace();
        }
        return null;
    }

生成流域tif文件的方法

 /**
     * 生成未裁剪的tif文件
     *
     * @param outTifFolder
     * @param shpFolder
     * @param ncDataset
     * @param geoTransform
     * @param nLat
     * @param nLon
     * @param driver
     * @param apcpArray
     * @param tifName
     * @throws IOException
     */
    private void genTif(String outTifFolder, String shpFolder, NetcdfFile ncDataset, double[] geoTransform,
                        int nLat, int nLon, Driver driver, float[][] apcpArray, String tifName) throws IOException {
        String fileAllPath = outTifFolder + GdalConstants.SPLIT + tifName;
        String fileAllPathTemp = outTifFolder + GdalConstants.SPLIT + "temp" + tifName;
        Dataset tifDataSet = driver.Create(fileAllPathTemp, nLon, nLat, 1, 6);

        tifDataSet.SetGeoTransform(geoTransform);
        //获取地理坐标系统信息, 定义输出的坐标系为"WGS 84"，AUTHORITY["EPSG","4326"]
        SpatialReference spatialReference = new SpatialReference();
        spatialReference.SetWellKnownGeogCS(GdalConstants.COORDINATE);
        String wkt = spatialReference.ExportToWkt();
        //给图层赋予投影信息
        tifDataSet.SetProjection(wkt);

        //更改异常值
        for (int i1 = 0; i1 < Objects.requireNonNull(apcpArray).length; i1++) {
            for (int i2 = 0; i2 < apcpArray[i1].length; i2++) {
                if (apcpArray[i1][i2] > 1000000) {
                    apcpArray[i1][i2] = 0f;
                }
            }
        }
        tifDataSet.GetRasterBand(1).SetNoDataValue(-32768);

        float[] arr = new float[nLon * nLat];
        int c = 0;
        for (int j = nLat - 1; j >= 0; j--) {
            for (int j2 = 0; j2 < nLon; j2++) {
                arr[c] = apcpArray[j][j2];
                c++;
            }
        }
        //参数说明： 起始x 起始y x数量(行) y数量(列) 数据类型 数组  //lat lon = 337 357 = x y
        tifDataSet.GetRasterBand(1).WriteRaster(0, 0, nLon, nLat, 6, arr);

        //重采样算法:三次卷积, 重采样参数设置
        Dataset tifDataSetReSam = driver.Create(fileAllPath, nLon * GdalConstants.RESAMPLE_FACTORY, nLat * GdalConstants.RESAMPLE_FACTORY, 1, 6);
        tifDataSetReSam.SetGeoTransform(getGeoTransform(ncDataset, GdalConstants.RESAMPLE_FACTORY));
        tifDataSetReSam.SetProjection(wkt);
        tifDataSetReSam.GetRasterBand(1).WriteRaster(0, 0, nLon, nLat, 6, arr);
        tifDataSetReSam.GetRasterBand(1).SetNoDataValue(GdalConstants.NO_DATA_VALUE);
        // 重采样
        gdal.ReprojectImage(tifDataSet, tifDataSetReSam, tifDataSet.GetProjection(), tifDataSetReSam.GetProjection(), 3, 0.0, 0.0);

        //写入硬盘
        tifDataSetReSam.FlushCache();
        //删除内存中的数据
        tifDataSet.delete();
        tifDataSetReSam.delete();
        //在生成tif之后, 需要按照指定的shape文件, 裁剪
        //按照需求暂时有两个shape需要裁剪, 后续会有接近20个shp
        List<File> shpPathList = new ArrayList<>();
        File[] files = new File(shpFolder).listFiles();
        if (Objects.isNull(files)) {
            throw new CustomException("当前shp文件夹为空");
        }
        for (File file : files) {
            if (file.getName().endsWith(".shp")) {
                shpPathList.add(file);
            }
        }
        if (shpPathList.size() == 0) {
            throw new CustomException("当前shp文件夹中没有shp文件,请检查");
        }
        tifClibByShp(fileAllPath, shpPathList, Boolean.TRUE);
        //删除没有重采样的数据
        new File(fileAllPathTemp).delete();
        //删除重采样后的大文件
        new File(fileAllPath).delete();
    }

置影像的显示范围

/**
     * 置影像的显示范围
     * -Lat_Res一定要是-的
     * 六个元素的元组:左上角x坐标，水平分辨率，旋转参数，左上角y坐标，旋转参数，竖直分辨率
     *
     * @param ncDataset      nc文件
     * @param resampleFactor 重采样因子
     * @return 影像的显示范围数组
     * @throws IOException
     */
    private static double[] getGeoTransform(NetcdfFile ncDataset, Integer resampleFactor) throws IOException {
        Variable latVar = ncDataset.findVariable(GdalConstants.LATITUDE);
        Variable lonVar = ncDataset.findVariable(GdalConstants.LONGITUDE);

        //把lon,lat,time数据处理成一个数组
        double[] latArray = (double[]) latVar.read().get1DJavaArray(Double.class);
        double[] lonArray = (double[]) lonVar.read().get1DJavaArray(Double.class);

        //影像的对角坐标,确定其位置
        Arrays.sort(latArray);
        double latMin = latArray[0];
        double latMax = latArray[latArray.length - 1];
        Arrays.sort(lonArray);
        double lonMin = lonArray[0];
        double lonMax = lonArray[lonArray.length - 1];

        //计算输出tif的分辨率
        //确定lat,lon都是一维数组
        int nLat = latArray.length;
        int nLon = lonArray.length;
        double lonRes = (lonMax - lonMin) / ((double) (nLon) - 1);
        double latRes = (latMax - latMin) / ((double) (nLat) - 1);
        //设置影像的显示范围 -Lat_Res一定要是-的 六个元素的元组:左上角x坐标，水平分辨率，旋转参数，左上角y坐标，旋转参数，竖直分辨率
        double[] geoTransform = {lonMin, lonRes / resampleFactor, 0, latMax, 0, -latRes / resampleFactor};
        return geoTransform;
    }

栅格裁剪的方法

private String tifClibByShp(String srcTif, List<File> shpPathList, Boolean isReSam) {


        File inputTif = new File(srcTif);
        String inputFolder = inputTif.getParent();
        String outTiffPath = "";

        for (File shpFile : shpPathList) {
            //LSH  TJH
            String shpName = shpFile.getName().split("\\.")[0];
            String inputTifName = shpName + "_" + inputTif.getName();
            //存储裁剪后的tif的文件夹
            String outTifFolder = inputFolder + GdalConstants.SPLIT + shpName;
            outTiffPath = outTifFolder + GdalConstants.SPLIT + inputTifName;

            File file = new File(outTifFolder);
            if (!file.exists()) {
                file.mkdir();
            }

            //使用执行命令的方式, 裁剪tif
            //gdalwarp -of GTiff -cutline LSH.shp 2022092508.003.tif src_dem.tif
            String command = String.format("gdalwarp -of GTiff -overwrite %s -cutline %s %s %s", isReSam ? "-crop_to_cutline" : "", shpFile.getAbsolutePath(), srcTif, outTiffPath);
            boolean flag = ProcessUtil.execCommand(command);
            if (flag) {
                log.info("{}, {} 裁剪成功", srcTif, shpFile.getName());
            } else {
                log.error("{}, {} 裁剪失败", srcTif, shpFile.getName());
            }
        }
        return outTiffPath;
    }