2013年1月12日土曜日

10mメッシュ標高データで 3D地図

プログラムは、「5mメッシュ標高データで3D地図」とほとんど同じです。
高さの比率 heightScale の値と、PerspectiveCamera の設定だけ変更しています。

高さの比率
 double heightScale = 20;

PerspectiveCameraの設定
<PerspectiveCamera 
                    FarPlaneDistance="1000" 
                    LookDirection="0,0,-1" 
                    UpDirection="0,1,0" 
                    NearPlaneDistance="2" 
                    Position="562,-375,700" 
                    FieldOfView="90" />

画像の場所は、「10mメッシュ標高データで地図を作成」のサンプルと同じで、沖縄県の辺野古周辺です。
FG-GML-3928-60-dem10b-20090201.xml


拡大図

5m メッシュ標高データで 3D地図

XMLデータのパースは、「5mメッシュ標高データで地図を作成」と同じです。
1つのメッシュをそれぞれ空間座標の1つの点としてみます。
ある点(x, y, z1)を左上とした4つの点からなる地形は、2つの三角形で構成します。
(x ,y ,z1)、(x + 1 ,y ,z2)、(x + 1 ,y + 1 ,z3) を頂点とする三角形と、(x ,y ,z1)、(x ,y + 1 ,z4)、(x + 1 ,y + 1 ,z3) を頂点とする三角形 です。
すべての点について、これらの三角形を、MeshGeometry3D の Positions に追加して、全体の地形を作っています。

「5mメッシュ標高データで地図を作成」のサンプルと同じ場所を3Dで表示してみました。
FG-GML-3927-15-02-DEM5B-20110915.xml


MainWindow.xaml

<Window x:Class="GisTest3.MainWindow"
        xmlns="http://schemas.microsoft.com/winfx/2006/xaml/presentation"
        xmlns:x="http://schemas.microsoft.com/winfx/2006/xaml"
        Height="350" Width="525">
    <Window.CommandBindings>
        <CommandBinding Command="Open" Executed="OpenCommandHandler"/>
    </Window.CommandBindings>
    <Grid>
        <Grid.RowDefinitions>
            <RowDefinition Height="Auto"/>
            <RowDefinition Height="*"/>
        </Grid.RowDefinitions>
        <Menu Grid.Row="0">
            <MenuItem Command="Open" />
        </Menu>
        <Viewport3D Grid.Row="1" Name="myViewport">
            <Viewport3D.Camera>
                <PerspectiveCamera 
                    FarPlaneDistance="200" 
                    LookDirection="0,0,-1" 
                    UpDirection="0,1,0" 
                    NearPlaneDistance="2" 
                    Position="112,-75,150" 
                    FieldOfView="90" />
            </Viewport3D.Camera>
            <Viewport3D.Children>
                <ModelVisual3D>
                    <ModelVisual3D.Content>
                        <Model3DGroup>
                            <Model3DGroup.Children>
                                <DirectionalLight Color="#FFFFFFFF" Direction="5,0,-4" />
                                <GeometryModel3D>
                                    <GeometryModel3D.Geometry>
                                        <MeshGeometry3D />
                                    </GeometryModel3D.Geometry>
                                    <GeometryModel3D.Material>
                                        <DiffuseMaterial>
                                            <DiffuseMaterial.Brush>
                                                <SolidColorBrush Color="White" />
                                            </DiffuseMaterial.Brush>
                                        </DiffuseMaterial>
                                    </GeometryModel3D.Material>
                                </GeometryModel3D>
                            </Model3DGroup.Children>
                        </Model3DGroup>
                    </ModelVisual3D.Content>
                </ModelVisual3D>
            </Viewport3D.Children>
        </Viewport3D>
    </Grid>
</Window>

MainWindow.xaml.cs

using Microsoft.Win32;
using System.Text.RegularExpressions;
using System.Windows;
using System.Windows.Input;
using System.Windows.Media.Media3D;
using System.Xml;

namespace GisTest3
{
    public partial class MainWindow Window
    {
        public MainWindow()
        {
            InitializeComponent();
        }

        private void OpenCommandHandler(object sender, ExecutedRoutedEventArgs e)
        {
            OpenFileDialog ofd = new OpenFileDialog();
            if (ofd.ShowDialog() == true)
            {
                XmlDocument xdoc = new XmlDocument();
                xdoc.Load(ofd.FileName);

                XmlNamespaceManager xnm = new XmlNamespaceManager(xdoc.NameTable);
                xnm.AddNamespace("gml""http://www.opengis.net/gml/3.2");

                XmlElement dataXe = (XmlElement)xdoc.SelectSingleNode("//gml:tupleList", xnm);
                XmlElement lowXe = (XmlElement)xdoc.SelectSingleNode("//gml:low", xnm);
                XmlElement highXe = (XmlElement)xdoc.SelectSingleNode("//gml:high", xnm);

                string[] lowPoint = lowXe.InnerText.Split(' ');
                string[] highPoint = highXe.InnerText.Split(' ');
                int xSize = int.Parse(highPoint[0]) - int.Parse(lowPoint[0]) + 1;
                int ySize = int.Parse(highPoint[1]) - int.Parse(lowPoint[1]) + 1;
                string data = dataXe.InnerText;

                DrawMap(xSize, ySize, data);
            }
        }

        private void DrawMap(int xSize, int ySize, string data)
        {
            float[,] heights = new float[xSize, ySize];

            Regex regex = new Regex("[^,\\r\\n]+,[^,\\r\\n]+");
            MatchCollection mc = regex.Matches(data);
            int index = 0;
            foreach (Match in mc)
            {
                string[] unit = m.Value.Split(',');
                heights[index % xSize, index / xSize] = (float.Parse(unit[1]) < 0) ? 0 : float.Parse(unit[1]);
                index++;
            }

            MeshGeometry3D meshGeometry3D = (MeshGeometry3D)((GeometryModel3D)((Model3DGroup)((ModelVisual3D)this.myViewport.Children[0]).Content).Children[1]).Geometry;
            meshGeometry3D.Positions.Clear();

            double heightScale = 3;

            for (int x = 0; x < xSize; x++)
            {
                for (int y = 0; y < ySize; y++)
                {
                    if (x < xSize - 1 && y < ySize - 1)
                    {
                        Point3D p11 = new Point3D(x, -1 * y, heights[x, y] / heightScale);
                        Point3D p12 = new Point3D(x + 1, -1 * y, heights[x + 1, y] / heightScale);
                        Point3D p13 = new Point3D(x + 1, -1 * (y + 1), heights[x + 1, y + 1] / heightScale);
                        Point3D p21 = new Point3D(x, -1 * y, heights[x, y] / heightScale);
                        Point3D p22 = new Point3D(x + 1, -1 * (y + 1), heights[x + 1, y + 1] / heightScale);
                        Point3D p23 = new Point3D(x, -1 * (y + 1), heights[x, y + 1] / heightScale);
                        meshGeometry3D.Positions.Add(p13);
                        meshGeometry3D.Positions.Add(p12);
                        meshGeometry3D.Positions.Add(p11);
                        meshGeometry3D.Positions.Add(p23);
                        meshGeometry3D.Positions.Add(p22);
                        meshGeometry3D.Positions.Add(p21);
                    }
                }
            }
        }
    }
}