Java读取气象水文Nc数据
最近业务上需要处理气象水文的NC数据，NC格式数据是气象水文存储的一种专业数据格式，也是GIS数据的一部分，今天我们就用Java来解析nc数据。
一、 环境准备
准备nc数据，搭建一个简单的maven环境即可。
二、解析nc
Java的三方库geotools的gt-netcdf库封装了对nc数据的解析代码，如果熟悉geotools朋友，可直接使用这个三方库，其实底层对nc数据的解析主要是基于如下的三方库
<dependency> <groupId>edu.ucar</groupId> <artifactId>cdm</artifactId> <version>4.6.15</version> </dependency> <dependency> <groupId>edu.ucar</groupId> <artifactId>netcdf4</artifactId> <version>4.6.15</version> </dependency> <dependency> <groupId>edu.ucar</groupId> <artifactId>bufr</artifactId> <version>4.6.15</version> </dependency>
因此，如果不想使用geotools，想自己解析，我们就需要将上面的依赖，添加到我们的pom中。具体如何解析呢，下面直接给出解析代码，然后会逐一说明nc的一些数据结构。
@Test public void renderChinaTemp() throws IOException, FactoryException, InvalidRangeException { //打开nc文件 NetcdfFile netcdfFile = NetcdfFile.open(FILE_PATH); //构造dataset NetcdfDataset netcdfDataset = new NetcdfDataset(netcdfFile); //坐标变量 List<CoordinateAxis> coordinateAxisList = netcdfDataset.getCoordinateAxes(); //获取nc数据的范围 Map<String,Object> map = getLonLatData(coordinateAxisList); //获取时间 getDataTimeList(coordinateAxisList,netcdfDataset); //请求的bbbox String bbox = "-180, 180, -90, 90";//获取其中的变量 List<Variable> variables = netcdfDataset.getVariables(); //获取无效值 double undefData = getUnvalidData(variables); /** * 获取格点数据,实际中第二个参数需要传递每个维度的下标 * 比如根据具体查询时间可以在上面的时间集合中获取到其对应的下标 * 对于经纬度的下标需要分别对待，一般纬度信息基本没有问题， * 范围基本都是在-90.90，我们只需要用查询的维度范围减去维度的起始值 * 的绝对值除以增幅就可以； * 对于经度范围一般全球都是0-360表示，如何将我们请求的-180—180，转到 * 0-360呢，其实呢只需要简单的转换下就行，只需要将<0的给他加360就行 * 也就是0-180不变对应0-360中的0-180 -180-0对应的是0——360中的180-360 * 其实个人理解：nc数据的采集时从本初子午线开始(0度经线，沿着自西向东的方向 * 先到+180°，后又从-180—0进行了采集) */ double[] latCoordinates = (double[]) map.get("latArray"); double[] lonCoordinates = (double[]) map.get("lonArray"); String resourceBbox = map.get("bbox").toString(); String[] bboxArray = resourceBbox.split(";"); double minLon = Double.parseDouble(bboxArray[0]); double maxLon = Double.parseDouble(bboxArray[1]); double minLat = Double.parseDouble(bboxArray[2]); double maxLat = Double.parseDouble(bboxArray[3]); //这里需要用原始请求范围和数据的范围求交集后在处理，由于我的数据和请求都是全球范围 Map<String,Object> paramMap = new HashMap<String,Object>(); //y方向 int starty,endy = -1; //这个主要是查看数据的记录是从高纬度开始，还是低纬度开始，用于计算下标 boolean upper = false; if (minLat > maxLat){ starty = (int) ((Math.abs(90 - minLat)) / 0.5); //其中0.5是数据的步幅，这个在计算范围的时候可以获取到 endy = (int) ((Math.abs(-90 - minLat)) / 0.5); }else{ starty = (int) ((Math.abs(-90 - minLat)) / 0.5); //其中0.5是数据的步幅，这个在计算范围的时候可以获取到 endy = (int) ((Math.abs(90 - minLat)) / 0.5); upper = true; } paramMap.put("lat",starty + ";" +(endy - 1)); //X方向，目前先处理请求范围0-180，根据上面的理论以及实际minLon-maxLon（0-360） int startx = (int) ((0 - minLon) / 0.5); int endx = (int) ((180 - minLon) / 0.5); paramMap.put("lon",startx + ";" + (endx - 1)); float[] precipArray = readNcData(variables,paramMap); }
nc数据主要参数是Variable变量，而这个变量有分为坐标维度变量和数据维度变量，一般我们拿到的nc数据，都是按天记录的，而且数据维度基本也是唯一的，但是如果我们从一些网站下载历史数据，可能会拿到时间维度不唯一，数据维度也不唯一。
坐标维度，一般数据中都是lon,lat,time三类，当然不同的数据维度，除了以上三个基本的维度外，还会有其它的维度信息，比如气压，高程等。
对于数据维度，这个也是核心，是一个nc数据存储数据的维度，比如温度，降雨量等，而对于这个维度，一般都会有其对应的Dimension，其实也就是上面说到的坐标维度，一般是通过坐标维度唯一确定一个格点数值，这样就可以通过一定的坐标维度，输出一个范围内的格点数据。
下面的代码是上面主类中的部分方法
private float getUnvalidData(List<Variable> variables){ for (Variable variable: variables) { //直接取其维度非坐标维度 if (!variable.isCoordinateVariable()) { //我们在这取第一个时间维度 List<Attribute> attributes = variable.getAttributes(); for (Attribute attribute : attributes) { String fullname = attribute.getFullName(); if ("missing_value".equalsIgnoreCase(fullname)) { return (float) attribute.getValue(0); } } } } return -99999.0f; } /** * 获取时间维度，并转换为Date * @param coordinateAxisList * @return */ public List<Date> getDataTimeList(List<CoordinateAxis> coordinateAxisList,NetcdfDataset netcdfDataset) throws IOException { for (CoordinateAxis coordinateAxis : coordinateAxisList) { AxisType axisType = coordinateAxis.getAxisType(); if (axisType.equals(AxisType.Time)) { List<Date> timeList = new ArrayList<>(); //获取时间维度，一般nc文件都是某一天的数据，但是有时候，会将数据合并,因此如果需要获取所有时间周期，需要处理如下逻辑 CoordinateAxis1DTime coordinateAxis1DTime = CoordinateAxis1DTime.factory(netcdfDataset, coordinateAxis, null); int[] shapeArray = coordinateAxis.getShape(); for (int i = 0; i < shapeArray[0]; i++) { Date date = coordinateAxis1DTime.getCalendarDate(i).toDate(); timeList.add(date); } return timeList; } } return null; } /** * 获取对应的bbox和经纬度变换范围 * @param coordinateAxisList * @return * @throws IOException */ public Map<String,Object> getLonLatData(List<CoordinateAxis> coordinateAxisList) throws IOException, FactoryException { float[] latArray = new float[0],lonArray = null; double[] latCoor = new double[2]; double[] lonCoor = new double[2]; Map<String,Object> result = new HashMap<String,Object>(); for (CoordinateAxis coordinateAxis : coordinateAxisList) { AxisType axisType = coordinateAxis.getAxisType(); //计算处坐标范围 if (axisType.equals(AxisType.Lat)){ getCoordinate((CoordinateAxis1D) coordinateAxis,latCoor); latArray = (float[]) coordinateAxis.read().copyTo1DJavaArray(); } if (axisType.equals(AxisType.Lon)){ getCoordinate((CoordinateAxis1D) coordinateAxis,lonCoor); lonArray = (float[]) coordinateAxis.read().copyTo1DJavaArray(); } } double[] doublelatArray = new double[latArray.length]; for (int i = 0; i < latArray.length; i++) { doublelatArray[i] = latArray[i]; } double[] doublelonArray = new double[lonArray.length]; for (int i = 0; i < lonArray.length; i++) { doublelonArray[i] = lonArray[i]; } result.put("latArray",doublelatArray); result.put("lonArray",doublelonArray); result.put("bbox",lonCoor[0]+";" + lonCoor[1] + ";" + latCoor[0] + ";" + latCoor[1]); return result; } private void getCoordinate(CoordinateAxis1D coordinateAxis1D,double[] coordinate){ if (coordinateAxis1D.isRegular()){ coordinate[0] = coordinateAxis1D.getStart() - (coordinateAxis1D.getIncrement() / 2d); coordinate[1] = coordinate[0] + coordinateAxis1D.getIncrement() * (coordinateAxis1D.getSize()); }else{ double min = ((Number) coordinateAxis1D.getMinValue()).doubleValue(); double max = ((Number) coordinateAxis1D.getMaxValue()).doubleValue(); double incr = (max - min) / (coordinateAxis1D.getSize() - 1); coordinate[0] = min - (incr / 2d); coordinate[1] = max + (incr / 2d); } } private float[] readNcData(List<Variable> variables, Map<String,Object> param) throws InvalidRangeException, IOException { float[] precipArray = null; for (Variable variable: variables) { //直接取其维度非坐标维度 if (!variable.isCoordinateVariable()){ //获取维度，查询每一个 List<Dimension> dimensions = variable.getDimensions(); Section section = new Section(); for (Dimension dimension : dimensions) { String dimensionName = dimension.getFullName(); if ("time".equalsIgnoreCase(dimensionName)) { //获取时间的下标维度 Range range = new Range(3, 3); section.appendRange(range); } if ("lon".equalsIgnoreCase(dimensionName)) { String lonRange = param.get("lon").toString(); String[] rangeArray = lonRange.split(";"); section.appendRange(Integer.parseInt(rangeArray[0]),Integer.parseInt(rangeArray[1])); } if ("lat".equalsIgnoreCase(dimensionName)){ //取全球数据，因此不做处理 String lonRange = param.get("lat").toString(); String[] rangeArray = lonRange.split(";"); section.appendRange(Integer.parseInt(rangeArray[0]),Integer.parseInt(rangeArray[1])); } } precipArray = (float[]) variable.read(section).copyTo1DJavaArray(); } } return precipArray; }
三、wContour生成等值线
上面已经解析了NC数据，并且生成了格点数据，对于业务而言，一般气象水文数据的指标，比如温度，降雨量，气压等都要生成等值线，因此在这里也顺便提及下nc数据如何利用wContour生成等值线，其实nc数据解析后已经是格点数据了，而且也有格点的分辨率，除非分辨率无法满足要求，否则nc的格点数据，就可以直接用来生成等值线。下面是生成等值线的代码。
/** * * @param _gridData nc的格点数据，y在前 x在后 * @param sampleInterval * @param _x 经度 * @param _y 维度 * @param undefData 无效数组 */ public void renderNcEquiLine(double[][] _gridData, double[] sampleInterval, double[] _x, double[] _y, double undefData) throws IOException { int[][] S1 = new int[_gridData.length][_gridData[0].length]; //获取轮廓边界 List<Border> _borders = Contour.tracingBorders(_gridData, _x, _y, S1, undefData); //生成等值线 List<PolyLine> cPolylineList = Contour.tracingContourLines(_gridData, _x, _y, sampleInterval.length, sampleInterval, undefData, _borders, S1); //将wcontour的polylin转化为geojson String lineGeojson = convertPolyline2GeoJson(cPolylineList); //将geojson转化为featurecolection SimpleFeatureCollection simpleFeatureCollection = convertGeojson2FeatureCollection(lineGeojson); //将SimpleFeatureCollection转化为shape ShapeFileWriter feShapeFileWriter = new ShapeFileWriter("D:\\Gis开发\\数据\\temp\\line.shp"); feShapeFileWriter.createScheam(simpleFeatureCollection.getSchema().getAttributeDescriptors(),"EPSG:4326"); feShapeFileWriter.writeShape(simpleFeatureCollection,simpleFeatureCollection.getSchema().getGeometryDescriptor().getLocalName(),true); //使用geotools的swing绘制查看效果 renderImg(simpleFeatureCollection); }
其中_gridData就是我们上面生成好的二维数据，但是记得这个二维数据是y在前，x在后的二维数据，并且数据的顺序是从左到右，从下向上的，因此如果NC存储的数据顺序与方法的不一致，需要进行反转，比如我的这个nc数据的lat是从上往下，因此调用这个等值线的方法时，需要将返回的一维数组转二维数组的时候，从下往上处理，具体代码如下：
//一维数组转为二维 double[][] gridData = convertOneArray2TwoArray(precipArray,(endy-starty),(endx - startx)); //等值线间隔 double[] sampInterval = {0,0.5,1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0,10,15}; double[] rangeCoord = Arrays.copyOfRange(lonCoordinates,0,360); double[] reverseCoord = new double[latCoordinates.length]; for (int i = 0;i<latCoordinates.length;i++){ reverseCoord[i] = latCoordinates[latCoordinates.length - 1 - i]; } //gridData数据一定要是从左倒有，从下到上的顺序 renderNcEquiLine(gridData,sampInterval,rangeCoord,reverseCoord,undefData);
好了，到这里nc数据的解析使用已经处理完毕了，如有问题，欢迎随时交流！对于wContour的使用，后面会专门出一篇wContour生成等值线和等值面的文章。
