java实现NC数据等值线等值面可视化

工作中的项目需要将NC科学数据进行等值线等值面处理，并展示给用户。网上找资料真费劲，基本没有直接可参考的，自己收集了些资料，参考着同事的思路实现效果了，跟大家分享下，也给自己做个记录。

刚开始写博客，写的不好的地方请大家见谅。

先上效果图。

由于数据边界比较敏感就给地图打马赛克了，截图也只有部分区域，不影响效果。

思路

1.首先java后端加载nc文件，读取nc变量的值，一般是多维数组。java读取NC文件用的是netcdf-4.3.22.jar。

2.通过第三方jar包对数组数据进行等值线处理,我用的wContour.jar,文末有下载链接。

3.后端生成的等值线数据返回到前端，通过leaflet展示。

java后端代码

主方法

	//生成等值线主方法
	public static Map<String, Object> isolineProcess() throws IOException{
		Map<String, Object> resMap = new HashMap<String, Object>();
		String filePath = "";//文件路径
		String element = "";//nc变量名
		int depthIndex = 1;//nc文件深度序列
		int timeIndex = 1;//nc文件时间序列
		//获取NC的数据
        Map map = getNcData(filePath,element,depthIndex,timeIndex);
        if(map == null || map.size()==0){
        	return resMap;
        }
        //等值线等值面间隔,间隔越小:过渡越均匀,性能消耗越大,响应越慢
        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, 35};
        String strGeojson = nc2EquiSurface(map, dataInterval);
        resMap.put("geojson", strGeojson);
        resMap.put("dataArr", map.get("eleData"));
        return resMap;
	}

处理nc：

// 获取nc数据
	public static <T> Map getNcData(String ncpath,String element,int depthIndex,int timeIndex) throws IOException {
		String lonVarName = "lon";//经度变量
        String latVarName = "lat";//纬度变量
		//加载nc文件
		NetcdfFile ncfile =  NetcdfDataset.open(ncpath);
		//读取经纬度数据
		Variable varLon = ncfile.findVariable(lonVarName);
		Variable varLat = ncfile.findVariable(latVarName);
		Object lonObj = null;
		Object latObj = null;
		lonObj = (Object) varLon.read().copyToNDJavaArray();
		latObj = (Object) varLat.read().copyToNDJavaArray();
		Map map = new HashMap();
		map = readNCLonLat(map,lonObj,latObj);
		int latLength = (int) map.get("latLength");
		int lonLength = (int) map.get("lonLength");
	
		//读取变量数据
		Variable varPre = ncfile.findVariable(element);
		int[] org = { timeIndex, depthIndex, 0, 0 };// 每一维的起点,起点从0开始
		int[] sha = { 1, 1, latLength, lonLength };// 相应维读取的条数
		/*
		 * nc数据类型可能是short、folat、double，所以用Object接收，通过反射判断数据类型
		 */
		Object pre = null;
		try {
			pre = (Object) varPre.read(org,sha).copyToNDJavaArray();
		} catch (InvalidRangeException e) {
			e.printStackTrace();
		}
		String className = pre.getClass().getComponentType().getName();
		
		Double invalidValueDou = -9999.0;//无效值
		Double scaleFactorDou = 1.0;//比例因子数值
		/*
		 * 这里只判断了short数组，其他数据类型的判断类似.
		 * 我的数据是四维的：经纬度，时间，深度。
		 */
		if("[[[s".equalsIgnoreCase(className)){
			short[][][][] dataArr = (short[][][][])pre;
			//创建新的double数组保存数据
			double[][] dPre = new double[dataArr[0][0].length][dataArr[0][0][0].length];
            //循环将short数组转为double数组
            for (int i = 0, len = dataArr[0][0].length; i < len; i++) {
            	short[] _pre = dataArr[0][0][i];
                for (int j = 0, jlen = _pre.length; j < jlen; j++) {
                	Double value = Double.parseDouble(String.valueOf(_pre[j]));
                	if(invalidValueDou != null && Double.compare(value,invalidValueDou)==0){
                		dPre[i][j] = invalidValueDou;
                		continue;//无效值跳过
                	}
                	if(scaleFactorDou == null){//比例因子处理
                		dPre[i][j] = value;
                	}else{
                		dPre[i][j] = value * scaleFactorDou;
                	}
                    
                }
            }
            map.put("eleData", dPre);
		}
		return map;
	}
	
	//处理nc经纬度
	private static Map readNCLonLat(Map map, Object lonObj, Object latObj) {
		/**
		 * nc数据类型可能是short、folat、double，所以用Object接收，通过反射判断数据类型,
		 * 这里只判断了short数组，其他数据类型的判断类似
		 */
		String lonclassName = lonObj.getClass().getComponentType().getName();
		int lonLength = 1;
		int latLength = 1;
		if ("[s".equalsIgnoreCase(lonclassName)) {
			short[] lonArr = (short[]) lonObj;
			short[] latArr = (short[]) latObj;
			lonLength = lonArr.length;
			latLength = latArr.length;
			//创建double数组存放经纬度数据，方便后续计算
			double[] dLon = new double[lonArr.length], dLat = new double[latArr.length];
			for (int i = 0, len = lonArr.length; i < len; i++) {
				dLon[i] = Double.parseDouble(String.valueOf(lonArr[i]));
			}
			for (int i = 0, len = latArr.length; i < len; i++) {
				dLat[i] = Double.parseDouble(String.valueOf(latArr[i]));
			}
			map.put("lon", dLon);
			map.put("lat", dLat);
		}
		return map;
	}

 

生成等值线：

//nc数据生成等值线数据
    public static String nc2EquiSurface(Map ncData, double[] dataInterval) {
        String geojsonpogylon = "";
 
        List<PolyLine> cPolylineList = new ArrayList<PolyLine>();
        List<Polygon> cPolygonList = new ArrayList<Polygon>();
 
        double[][] _gridData = (double[][]) ncData.get("eleData");
        int[][] S1 = new int[_gridData.length][_gridData[0].length];
        double[] _X = (double[]) ncData.get("lon"), _Y = (double[]) ncData.get("lat");
        double _undefData = Double.parseDouble((String)ncData.get("invalid"));
        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);// 生成等值线
 
        cPolylineList = Contour.smoothLines(cPolylineList);// 平滑
        cPolygonList = Contour.tracingPolygons(_gridData, cPolylineList,
                _borders, dataInterval);
 
        geojsonpogylon = getPolygonGeoJson(cPolygonList);
 
        return geojsonpogylon;
    }
    
    //拼装geogson
    public static String getPolygonGeoJson(List<Polygon> cPolygonList) {
		String geo = null;
		String geometry = " { \"type\":\"Feature\",\"geometry\":";
		String properties = ",\"properties\":{ \"value\":";
		
		String head = "{\"type\": \"FeatureCollection\"," + "\"features\": [";
		String end = "  ] }";
		if (cPolygonList == null || cPolygonList.size() == 0) {
			return null;
		}
		try {
			for (Polygon pPolygon : cPolygonList) {
				List<Object> ptsTotal = new ArrayList<Object>();
					for (PointD ptd : pPolygon.OutLine.PointList) {
						List<Double> pt = new ArrayList<Double>();
						pt.add(doubleFormat(ptd.X));
	                    pt.add(doubleFormat(ptd.Y));
						ptsTotal.add(pt);
					}
				List<Object> list3D = new ArrayList<Object>();
				list3D.add(ptsTotal);
				JSONObject js = new JSONObject();
				js.put("type", "Polygon");
				js.put("coordinates", list3D);
				
				geo = geometry + js.toString() + properties  +pPolygon.LowValue + "} }" + "," + geo;
			}
			if (geo.contains(",")) {
				geo = geo.substring(0, geo.lastIndexOf(","));
			}
			
			geo = head + geo + end;
		} catch (Exception e) {
			e.printStackTrace();
			return geo;
		}
		return geo;
	}
    
    /**
   	 * double保留两位小数
   	 */
   	public static double doubleFormat(double d) {  
           BigDecimal bg = new BigDecimal(d);  
           double f1 = bg.setScale(4, BigDecimal.ROUND_HALF_UP).doubleValue();  
           return f1;
       } 

后端代码基本就这些，最后返回的是geoJson格式的数据。后端处理过程还是比较慢的，我这个效果本地测试大约3秒，如果用MATLAB做应该会好些，但是不会MATLAB~~~

其实到这基本就可以了，前端拿着数据渲染就行，我是用leaflet实现的。

等之后有时间了在把前端js代码补充下。

另外还有一种方法是把nc数据直接返回到前端，js进行等值线处理，这种方法就是比较消耗浏览器性能。

 

---

2020-12-12更新,

netcdf-4.3.22.jar和wContour.jar上传了,需要的朋友点击下载

把前端代码补上: 前端通过ajax获取到数据后渲染

// 添加等值线,等值面
function loadLine(geojsonData) {
	
	var maxTemp = 35;
	var minTemp = 0;
	var rgbValue = "";
	$.getJSON('data/json/color.json',function(colorJson){
		let num = colorJson.colorList.length;
		arr = geojsonData.features;
		// 循环添加多边形
		for (var i = 0; i < arr.length; i++) {
			var dataone = arr[i].geometry.coordinates[0];
			for (var j = 0; j < dataone.length; j++) {
				var temp1 = dataone[j][0];
				var temp2 = dataone[j][1];
				dataone[j][0] = temp2;
				dataone[j][1] = temp1;
			}
			//var value = Math.round(arr[i].properties.value);
			var value = arr[i].properties.value;
			// console.log(value);
			//获取颜色rgb.颜色赋值
			var rgb = "";
			if(value >= 30){  //超过30°取最后一个
				rgb =  "rgb("+colorJson.colorList[212].color213+")";
			}else if(value <= 0){  //低于0°取第一个
				rgb =  "rgb("+colorJson.colorList[0].color1+")";
			}else{  //0°~30°
			let interval = 30/num; //最高30°，最低0°
			let index = parseInt(value/interval);
			let colorObj = colorJson.colorList[index-1];
			let name = "color"+index;
			rgb =  "rgb("+colorObj[name]+")";
			}
			polygon = L.polygon(dataone).addTo(polygon_group);
			polygon.setStyle({
				color : rgb,// 
				//color :"rgb(0,192,255)",
				weight : 0,
				fillOpacity : 1
			});
		}
		
	});
	//创建图层,设置层级
	map.createPane('ncline');
	map.getPane('ncline').style.zIndex = 800;
	L.geoJson(geojsonData, {
		 color: "black",
		 fill:false,
		 pane: 'ncline',
		 weight: 1,
		 lineCap:"5",
		 text:"25"
	}).addTo(polygon_group);
	polygon_group.addTo(map);
}

 

总结:上述代码等值线效果还行,不过最近用了geoserver,感觉展示效果也不错,大家也可以试试