1 <dependency>
2 <groupId>edu.ucar</groupId>
3 <artifactId>netcdfAll</artifactId>
4 <version>5.5.4</version>
5 <scope>system</scope>
6 <systemPath>${pom.basedir}/../lib/netcdfAll-5.5.4-SNAPSHOT-x.jar</systemPath>
7 </dependency>
8 <dependency>
9 <groupId>org.meteothink</groupId>
10 <artifactId>wContour</artifactId>
11 <version>1.7.1</version>
12 </dependency>
1 2
3 import org.json.JSONObject;
4 import ucar.nc2.NetcdfFile;
5 import ucar.nc2.Variable;
6 import ucar.nc2.dataset.NetcdfDataset;
7 import wcontour.Contour;
8 import wcontour.global.Border;
9 import wcontour.global.PointD;
10 import wcontour.global.PolyLine;
11 import wcontour.global.Polygon;
12
13 import java.io.BufferedWriter;
14 import java.io.FileWriter;
15 import java.io.IOException;
16 import java.math.BigDecimal;
17 import java.text.DecimalFormat;
18 import java.util.*;
19
20 public class IsosurfacesUtil {
21
22 private static String rootPath = "E:\\Desktop\\";
23
24 private static final double _undefData = -9999.0;
25 private static final double invalidValueDou = -9999.0;
26 private static final double scaleFactorDou = 1.0;
27
28 public static void main(String[] args) throws IOException {
29 Map<String, Object> stringObjectMap = isolineProcess();
30 String strGeojson = (String) stringObjectMap.get("geojson");
31 String strFile = rootPath + "012.geojson";
32 BufferedWriter writer = new BufferedWriter(new FileWriter(strFile));
33 writer.write(strGeojson);
34 writer.close();
35 }
36
37
38
39
40 public static Map<String, Object> isolineProcess() throws IOException {
41 Map<String, Object> resMap = new HashMap<String, Object>();
42 String filePath = "E:\\Desktop\\012.grib";
43 String element = "Temperature_surface";
44 int depthIndex = 2;
45 int timeIndex = 1;
46 //获取NC的数据
47 Map map = getNcData(filePath,element,depthIndex,timeIndex);
48 if(map == null || map.size()==0){
49 return resMap;
50 }
51 double[] dataInterval = new double[]{-14,-13,-12,-11,-10,-9,-8,-7,-6,-5,-4,-3,-2,-1,0,1,2,3,4,5,6,7,8,9,10,11,12,13,14};
52 String strGeojson = nc2EquiSurface(map, dataInterval);
53 resMap.put("geojson", strGeojson);
54 resMap.put("dataArr", map.get("eleData"));
55 return resMap;
56 }
57
58
59 public static <T> Map getNcData(String ncpath,String element,int depthIndex,int timeIndex) throws IOException {
60 String lonVarName = "lon";
61 String latVarName = "lat";
62 //加载nc文件
63 NetcdfFile ncfile = NetcdfDataset.open(ncpath);
64 //读取经纬度数据
65 Variable varLon = ncfile.findVariable(lonVarName);
66 Variable varLat = ncfile.findVariable(latVarName);
67 Object lonObj = null;
68 Object latObj = null;
69 lonObj = (Object) varLon.read().copyToNDJavaArray();
70 latObj = (Object) varLat.read().copyToNDJavaArray();
71 Map map = new HashMap();
72 map = readNCLonLat(map,lonObj,latObj);
73 int latLength = (int) map.get("latlength");
74 int lonLength = (int) map.get("lonlength");
75
76 //读取变量数据
77 Variable varPre = ncfile.findVariable(element);
78
79 Object pre = null;
80 pre = (Object) varPre.read().copyToNDJavaArray();
81 String className = pre.getClass().getComponentType().getName();
82
83 if("[[f".equalsIgnoreCase(className)){
84 float[][][] dataArr = (float[][][])pre;
85 double[][] dPre = new double[dataArr[0].length][dataArr[0][0].length];
86 DecimalFormat df = new DecimalFormat("#.00");
87 for (int i = 0, len = dataArr[0].length; i < len; i++) {
88 float[] _pre = dataArr[0][i];
89 for (int j = _pre.length - 1, jlen = 0; j >= jlen; j--) {
90 double value = Double.parseDouble(String.valueOf(_pre[j]))- 273.15;
91 if(Double.compare(value, invalidValueDou) == 0){
92 dPre[i][j] = invalidValueDou;
93 continue;//无效值跳过
94 }
95 dPre[i][j] = Double.parseDouble(df.format(value * scaleFactorDou));
96
97 }
98 }
99 map.put("eleData", dPre);
100 }
101 return map;
102 }
103
104
105 private static Map readNCLonLat(Map map, Object lonObj, Object latObj) {
106
107 String lonclassName = lonObj.getClass().getComponentType().getName();
108 int lonLength = 1;
109 int latLength = 1;
110 if ("[s".equalsIgnoreCase(lonclassName)) {
111 short[] lonArr = (short[]) lonObj;
112 short[] latArr = (short[]) latObj;
113 lonLength = lonArr.length;
114 latLength = latArr.length;
115 //创建double数组存放经纬度数据，方便后续计算
116 double[] dLon = new double[lonArr.length], dLat = new double[latArr.length];
117 for (int i = 0, len = lonArr.length; i < len; i++) {
118 dLon[i] = Double.parseDouble(String.valueOf(lonArr[i]));
119 }
120 for (int i = 0, len = latArr.length - 1 ; i >= len; i--) {
121 dLat[i] = Double.parseDouble(String.valueOf(latArr[latArr.length-i-1]));
122 }
123 map.put("lon", dLon);
124 map.put("lat", dLat);
125 }
126 if ("float".equalsIgnoreCase(lonclassName)) {
127 float[] lonArr = (float[]) lonObj;
128 float[] latArr = (float[]) latObj;
129 lonLength = lonArr.length;
130 latLength = latArr.length;
131 //创建double数组存放经纬度数据，方便后续计算
132 double[] dLon = new double[lonArr.length], dLat = new double[latArr.length];
133 for (int i = 0, len = lonArr.length; i < len; i++) {
134 dLon[i] = Double.parseDouble(String.valueOf(lonArr[i]));
135 }
136 for (int i = 0, len = latArr.length; i < len; i++) {
137 dLat[i] = Double.parseDouble(String.valueOf(latArr[latArr.length-i-1]));
138 }
139 map.put("lon", dLon);
140 map.put("lat", dLat);
141 }
142 map.put("lonlength", lonLength);
143 map.put("latlength", latLength);
144 return map;
145 }
146
147
148 public static String nc2EquiSurface(Map ncData, double[] dataInterval) {
149 String geojsonpogylon = "";
150 String geojsonpogylona = "";
151
152 List<PolyLine> cPolylineList;
153 List<Polygon> cPolygonList;
154
155 double[][] _gridData = (double[][]) ncData.get("eleData");
156 int[][] S1 = new int[_gridData.length][_gridData[0].length];
157 double[] _X = (double[]) ncData.get("lon"), _Y = (double[]) ncData.get("lat");
158 System.out.println("dataInterval: " + Arrays.toString(dataInterval));
159 System.out.println("Min temperature: " + Arrays.stream(_gridData).flatMapToDouble(Arrays::stream).min().getAsDouble());
160 System.out.println("Max temperature: " + Arrays.stream(_gridData).flatMapToDouble(Arrays::stream).max().getAsDouble());
161 System.out.println("_gridData dimensions: " + _gridData.length + " x " + _gridData[0].length);
162 System.out.println("_X length: " + _X.length);
163 System.out.println("_Y length: " + _Y.length);
164
165 List<Border> _borders = Contour.tracingBorders(_gridData, _X, _Y, S1, _undefData);
166 int nc = dataInterval.length;
167 cPolylineList = Contour.tracingContourLines(_gridData, _X, _Y, nc, dataInterval, _undefData, _borders, S1);// 生成等值线
168
169 cPolylineList = Contour.smoothLines(cPolylineList);// 平滑
170 cPolygonList = Contour.tracingPolygons(_gridData, cPolylineList,_borders, dataInterval);
171
172 geojsonpogylon = getPolygonGeoJson(cPolygonList);
173
174 // geojsonpogylona = getPolylineGeoJson(cPolylineList);
175 return geojsonpogylon;
176 }
177
178
179
180 public static String getPolygonGeoJson(List<Polygon> cPolygonList) {
181 String geo = null;
182 String geometry = " { \"type\":\"Feature\",\"geometry\":";
183 String properties = ",\"properties\":{ \"value\":";
184
185 String head = "{\"type\": \"FeatureCollection\"," + "\"features\": [";
186 String end = " ] }";
187 if (cPolygonList == null || cPolygonList.size() == 0) {
188 return null;
189 }
190 try {
191 for (Polygon pPolygon : cPolygonList) {
192 List<Object> ptsTotal = new ArrayList<Object>();
193 for (PointD ptd : pPolygon.OutLine.PointList) {
194 List<Double> pt = new ArrayList<Double>();
195 pt.add(doubleFormat(ptd.X));
196 pt.add(doubleFormat(ptd.Y));
197 ptsTotal.add(pt);
198 }
199 List<Object> list3D = new ArrayList<Object>();
200 list3D.add(ptsTotal);
201 JSONObject js = new JSONObject();
202 js.put("type", "Polygon");
203 js.put("coordinates", list3D);
204
205 geo = geometry + js + properties +pPolygon.LowValue + "} }" + "," + geo;
206 }
207 if (geo.contains(",")) {
208 geo = geo.substring(0, geo.lastIndexOf(","));
209 }
210
211 geo = head + geo + end;
212 } catch (Exception e) {
213 e.printStackTrace();
214 return geo;
215 }
216 return geo;
217 }
218
219 public static String getPolylineGeoJson(List<PolyLine> cPolylineList) {
220 if (cPolylineList == null || cPolylineList.size() == 0) {
221 return null;
222 }
223
224 String collectionFormat = "{\"type\": \"FeatureCollection\"," + "\"features\": [ %s ] }";
225 String geometryFormat = "{\"type\":\"Feature\",\"geometry\": %s ,\"properties\":{ \"value\": %s }}";
226
227 List<String> geoList = new ArrayList<>();
228 for (int i = 0; i < cPolylineList.size(); i++) {
229 PolyLine pPolyline = cPolylineList.get(i);
230 List<Object> ptsTotal = new ArrayList<Object>();
231 for (PointD ptd : pPolyline.PointList) {
232 List<Double> pt = new ArrayList<Double>();
233 pt.add(doubleFormat(ptd.X));
234 pt.add(doubleFormat(ptd.Y));
235 ptsTotal.add(pt);
236 }
237 cn.hutool.json.JSONObject js = new cn.hutool.json.JSONObject();
238 js.set("type", "LineString");
239 js.set("coordinates", ptsTotal);
240 js.set("id", pPolyline.Value + "_" + i);
241
242 String geo = String.format(geometryFormat, js, pPolyline.Value);
243 geoList.add(geo);
244 }
245 return String.format(collectionFormat, String.join(", ", geoList));
246 }
247
248 /**
249 * double保留两位小数
250 */
251 public static double doubleFormat(double d) {
252 BigDecimal bg = new BigDecimal(d);
253 double f1 = bg.setScale(4, BigDecimal.ROUND_HALF_UP).doubleValue();
254 return f1;
255 }
256 }
浙公网安备 33010602011771号