java读取nc生成等值线

依赖

		<dependency>
			<groupId>edu.ucar</groupId>
			<artifactId>netcdf</artifactId>
			<version>4.2.20</version>
		</dependency>
		<dependency>
			<groupId>org.meteothink</groupId>
			<artifactId>wContour</artifactId>
			<version>1.7.1</version>
		</dependency>
		<dependency>
			<groupId>cn.hutool</groupId>
			<artifactId>hutool-all</artifactId>
			<version>5.8.32</version>
		</dependency>

代码：

package com.test.nc;

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.math.BigDecimal;
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

import cn.hutool.json.JSONObject;
import ucar.nc2.NetcdfFile;
import ucar.nc2.Variable;
import ucar.nc2.constants.AxisType;
import ucar.nc2.dataset.CoordinateAxis;
import ucar.nc2.dataset.NetcdfDataset;
import wcontour.Contour;
import wcontour.global.Border;
import wcontour.global.PointD;
import wcontour.global.PolyLine;
import wcontour.global.Polygon;

public class Nc2EquiSurface {
	private static String rootPath = System.getProperty("user.dir");
	private static double _undefData = -99999.0;
	
	public Map getNcData(String ncpath) throws Exception {
		Map map = new HashMap();
		NetcdfFile ncfile = NetcdfFile.open(ncpath);
		// 注意变量名字和类型，都转成double，因为wcontour需要
		Variable varLon = ncfile.findVariable("lon");
		Variable varLat = ncfile.findVariable("lat");
		Variable varPre = ncfile.findVariable("z");
		double[] dLon = (double[]) varLon.read().copyToNDJavaArray();
		double[] dLat = (double[]) varLat.read().copyToNDJavaArray();
		float[][] pre = (float[][]) varPre.read().copyToNDJavaArray();
		double[][] dPre = new double[pre.length][pre[0].length];
		for (int i = 0, len = pre.length; i < len; i++) {
			float[] _pre = pre[i];
			for (int j = 0, jlen = _pre.length; j < jlen; j++) {
				if (String.valueOf(_pre[j]).equals("NaN"))
					dPre[i][j] = _undefData;
				else
					dPre[i][j] = Double.parseDouble(String.valueOf(_pre[j]));
			}
		}

		map.put("lon", dLon);
		map.put("lat", dLat);
		map.put("tem", dPre);
		return map;
	}

	public String nc2EquiSurface(Map ncData, double[] dataInterval) {

		List<PolyLine> cPolylineList = new ArrayList<PolyLine>();
		List<Polygon> cPolygonList = new ArrayList<Polygon>();

		double[][] _gridData = (double[][]) ncData.get("tem");
		int[][] S1 = new int[_gridData.length][_gridData[0].length];
		double[] _X = (double[]) ncData.get("lon"), _Y = (double[]) ncData.get("lat");

		List<Border> _borders = Contour.tracingBorders(_gridData, _X, _Y, S1, _undefData);

		int nc = dataInterval.length;
		cPolylineList = Contour.tracingContourLines(_gridData, _X, _Y, nc, dataInterval, _undefData, _borders, S1);// 生成等值线
		return getPolylineGeoJson(cPolylineList);
	}
	
	public static String getPolylineGeoJson(List<PolyLine> cPolylineList) {
		if (cPolylineList == null || cPolylineList.size() == 0) {
			return null;
		}

		String collectionFormat = "{\"type\": \"FeatureCollection\"," + "\"features\": [ %s ] }";
		String geometryFormat = "{\"type\":\"Feature\",\"geometry\": %s ,\"properties\":{ \"value\": %s }}";

		List<String> geoList = new ArrayList<>();
		for (int i = 0; i < cPolylineList.size(); i++) {
			PolyLine pPolyline = cPolylineList.get(i);
			List<Object> ptsTotal = new ArrayList<Object>();
			for (PointD ptd : pPolyline.PointList) {
				List<Double> pt = new ArrayList<Double>();
				pt.add(doubleFormat(ptd.X));
				pt.add(doubleFormat(ptd.Y));
				ptsTotal.add(pt);
			}
			JSONObject js = new JSONObject();
			js.set("type", "LineString");
			js.set("coordinates", ptsTotal);
			js.set("id", pPolyline.Value + "_" + i); // 设置id

			String geo = String.format(geometryFormat, js.toString(), pPolyline.Value);
			geoList.add(geo); // 添加到结尾
		}
		return String.format(collectionFormat, String.join(", ", geoList));
	}
	
    // 截取坐标长度，过长没啥用
	public static double doubleFormat(double d) {
		BigDecimal bg = new BigDecimal(d);
		double f1 = bg.setScale(4, BigDecimal.ROUND_HALF_UP).doubleValue();
		return f1;
	}

	public static void main(String[] args) throws Exception {
		Nc2EquiSurface nc2equ = new Nc2EquiSurface();

		long start = System.currentTimeMillis();
		String ncpath = rootPath + "/data/ydy.nc";

		// 获取NC的数据
		Map map = nc2equ.getNcData(ncpath);
		
		double[] dataInterval = new double[] { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,
				18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30 };
		String strGeojson = nc2equ.nc2EquiSurface(map, dataInterval);

		String strFile = rootPath + "/data/ydy.geojson";
		BufferedWriter writer = new BufferedWriter(new FileWriter(strFile));
		writer.write(strGeojson);
		writer.close();

		System.out.println("Total cost:" + (System.currentTimeMillis() - start));
	}
}

后续需要优化：

1、生成的等值线数据量大，利用geotools简化。

2、wcontour批量生成的可能会缺数据，改成一条一条生成再合并。