WEBGIS开发教程

Appearance

章节导航

等值面生成 ​

一、wcontour简介 ​

wContour类库是使用 Java 和 C# 开发的，用于在使用 Java 或 .NET 环境的地球科学数据分析软件开发中使用等高线相关算法。对具有大量不规则分布的未定义值点的复杂网格数据的等高线追踪、填充和裁剪算法进行了全面的描述。从头开始描绘边界线，最后输出带孔的轮廓多边形。为网格数据的离散数据插值提供了反距离加权 (IDW) 和 Cressman 客观分析方法。包括轮廓线跟踪、平滑、填充和裁剪算法，以实现复杂网格数据的轮廓功能。还提供了流线分析方法。

二、实现 ​

1. 训练数据 ​

数据格式为csv的全国的3000多个气象站的温度数据，格式如下：

121.474042,30.917795,22
109.879153,31.074834,18
106.551541,29.563121,27
......

2. 数据解析与处理 ​

解析文件，并生成trainData。

java
String csvpath = "./tem.csv";
File csv = new File(csvpath);
BufferedReader br = null;
try {
	List<Double> lon = new ArrayList<Double>(),
         lat = new ArrayList<Double>(),
         heat = new ArrayList<Double>();
	br = new BufferedReader(new FileReader(csv));
	String line = "";
	while ((line = br.readLine()) != null) {
		String[] lineData = line.split(",");
		lon.add(Double.parseDouble(lineData[0]));
		lat.add(Double.parseDouble(lineData[1]));
		heat.add(Double.parseDouble(lineData[2]));
	}

	double[][] trainData = new double[3][lon.size()];
	for (int i = 0; i < lon.size(); i++) {
		trainData[0][i] = lon.get(i);
		trainData[1][i] = lat.get(i);
		trainData[2][i] = heat.get(i);
	}

} catch (Exception e) {
  e.printStackTrace();
}

3. 插值计算 ​

输入三个参数：

  • trainData，训练数据
  • dataInterval，数据间隔
java
double[] dataInterval = new double[]{5, 10, 15, 20, 25, 30, 35, 40, 45}
  • size，网格大小
java
int[] size = new int[]{100, 100};

如下图： 等值面

  • boundryFile，边界shp文件，一方面用于计算插值的范围，一方面如有需要，用于裁剪插值结果。 插值计算的代码如下：
java
import wContour.Contour;
import wContour.Global.Border;
import wContour.Global.PointD;
import wContour.Global.PolyLine;
import wContour.Global.Polygon;
import wContour.Interpolate;

try {
    double _undefData = -9999.0;
    SimpleFeatureCollection polygonCollection = null;
    List<PolyLine> cPolylineList = new ArrayList<PolyLine>();
    List<Polygon> cPolygonList = new ArrayList<Polygon>();

    int width = size[0],
            height = size[1];
    double[] _X = new double[width];
    double[] _Y = new double[height];

    File file = new File(boundryFile);
    ShapefileDataStore shpDataStore = null;

    shpDataStore = new ShapefileDataStore(file.toURL());
    //设置编码
    Charset charset = Charset.forName("GBK");
    shpDataStore.setCharset(charset);
    String typeName = shpDataStore.getTypeNames()[0];
    SimpleFeatureSource featureSource = null;
    featureSource = shpDataStore.getFeatureSource(typeName);
    SimpleFeatureCollection fc = featureSource.getFeatures();

    double minX = fc.getBounds().getMinX();
    double minY = fc.getBounds().getMinY();
    double maxX = fc.getBounds().getMaxX();
    double maxY = fc.getBounds().getMaxY();

    Interpolate.CreateGridXY_Num(minX, minY, maxX, maxY, _X, _Y);
    double[][] _gridData = new double[width][height];

    int nc = dataInterval.length;

    _gridData = Interpolate.Interpolation_IDW_Neighbor(trainData,
            _X, _Y, 12, _undefData);// IDW插值

    int[][] S1 = new int[_gridData.length][_gridData[0].length];
    /**
     * double[][] S0,
     * double[] X,
     * double[] Y,
     * int[][] S1,
     * double undefData
     */
    List<Border> _borders = Contour.tracingBorders(_gridData, _X, _Y,
            S1, _undefData);

    /**
     * double[][] S0,
     * double[] X,
     * double[] Y,
     * int nc,
     * double[] contour,
     * double undefData,
     * List<Border> borders,
     * int[][] S1
     */
    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);

    if (isclip) {
        polygonCollection = GeoJSONUtil.readGeoJsonByString(geojsonpogylon);
        FeatureSource dc = clipFeatureCollection(fc, polygonCollection);
        geojsonpogylon = getPolygonGeoJson(dc.getFeatures());
    }
} catch (Exception e) {
    e.printStackTrace();
}

4. 裁剪 ​

java
public FeatureSource clipFeatureCollection(FeatureCollection fc, SimpleFeatureCollection gs) {
    FeatureSource cs = null;
    try {
        List<Map<String, Object>> values = new ArrayList<Map<String, Object>>();
        FeatureIterator contourFeatureIterator = gs.features();
        FeatureIterator dataFeatureIterator = fc.features();
        while (dataFeatureIterator.hasNext()) {
            Feature dataFeature = dataFeatureIterator.next();
            Geometry dataGeometry = (Geometry) dataFeature.getProperty(
                    "the_geom").getValue();
            while (contourFeatureIterator.hasNext()) {
                Feature contourFeature = contourFeatureIterator.next();
                Geometry contourGeometry = (Geometry) contourFeature
                        .getProperty("geometry").getValue();
                double lv = (Double) contourFeature.getProperty("lvalue")
                        .getValue();
                double hv = (Double) contourFeature.getProperty("hvalue")
                        .getValue();
                if (dataGeometry.intersects(contourGeometry)) {
                    Geometry geo = dataGeometry
                            .intersection(contourGeometry);
                    Map<String, Object> map = new HashMap<String, Object>();
                    map.put("the_geom", geo);
                    map.put("lvalue", lv);
                    map.put("hvalue", hv);
                    values.add(map);
                }

            }

        }

        contourFeatureIterator.close();
        dataFeatureIterator.close();

        SimpleFeatureCollection sc = FeaureUtil
                .creatSimpleFeatureByFeilds(
                        "polygons",
                        "crs:4326,the_geom:MultiPolygon,lvalue:double,hvalue:double",
                        values);
        cs = FeaureUtil.creatFeatureSourceByCollection(sc);

    } catch (Exception e) {
        e.printStackTrace();
        return cs;
    }

    return cs;
}

生成后的结果如下如： 等值面

Copyright © 2023 LZUGIS