Skip to content
Navigation Menu
{{ message }}
forked from SharpMap/SharpMap
-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathProjectionExamples.cs
More file actions
472 lines (384 loc) · 19.1 KB
/
Copy pathProjectionExamples.cs
File metadata and controls
472 lines (384 loc) · 19.1 KB
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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
//#define alglib
namespace ExampleCodeSnippets
{
[NUnit.Framework.TestFixture]
public class ProjectionExamples
{
private const string Osgb36 =
"COMPD_CS[\"OSGB36 / British National Grid + ODN\",PROJCS[\"OSGB 1936 / British National Grid\",GEOGCS[\"OSGB 1936\",DATUM[\"OSGB 1936\",SPHEROID[\"Airy 1830\",6377563.396,299.3249646,AUTHORITY[\"EPSG\",\"7001\"]],TOWGS84[446.448,-125.157,542.06,0.15,0.247,0.842,-4.2261596151967575],AUTHORITY[\"EPSG\",\"6277\"]],PRIMEM[\"Greenwich\",0.0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.017453292519943295],AXIS[\"Geodetic latitude\",NORTH],AXIS[\"Geodetic longitude\",EAST],AUTHORITY[\"EPSG\",\"4277\"]],PROJECTION[\"Transverse Mercator\",AUTHORITY[\"EPSG\",\"9807\"]],PARAMETER[\"central_meridian\",-2.0],PARAMETER[\"latitude_of_origin\",49.0],PARAMETER[\"scale_factor\",0.9996012717],PARAMETER[\"false_easting\",400000.0],PARAMETER[\"false_northing\",-100000.0],UNIT[\"m\",1.0],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"27700\"]],VERT_CS[\"Newlyn\",VERT_DATUM[\"Ordnance Datum Newlyn\",2005,AUTHORITY[\"EPSG\",\"5101\"]],UNIT[\"m\",1.0],AXIS[\"Gravity-related height\",UP],AUTHORITY[\"EPSG\",\"5701\"]],AUTHORITY[\"EPSG\",\"7405\"]]";
private const string WGS84 =
"GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.01745329251994328,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]]";
[NUnit.Framework.Test]
[NUnit.Framework.Ignore("Nothing to see here")]
public void TestConversionProjNet()
{
var csf = new ProjNet.CoordinateSystems.CoordinateSystemFactory();
var cs1 = csf.CreateFromWkt(Osgb36);
var cs2 = csf.CreateFromWkt(WGS84);
var ctf = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();
var ct = ctf.CreateFromCoordinateSystems(cs1, cs2);
System.Diagnostics.Debug.Assert(ct != null);
}
#if alglib
/// <summary>
/// Performs an affine 2D coordinate transfromation
/// X' = _a*X + _b*Y + _c
/// Y' = _d*X + _e*Y + _f
/// </summary>
public class AffineCoordinateTransformation2D : ProjNet.CoordinateSystems.Transformations.MathTransform, GeoAPI.CoordinateSystems.Transformations.ICoordinateTransformation
{
private readonly double _a, _b, _c, _d, _e, _f;
private readonly double _ainv, _binv, _cinv, _dinv, _einv, _finv;
//private Matrix matrix;
//private readonly Matrix _inverse;
/// <summary>
/// Creates an instance of this class
/// </summary>
/// <param name="matrix"><see cref="T:System.Drawing.Drawing2D.Matrix"/> that holds the coefficients and for the coordinate transformation</param>
public AffineCoordinateTransformation2D(System.Drawing.Drawing2D.Matrix matrix)
: this(
matrix.Elements[0], matrix.Elements[1], matrix.OffsetX,
matrix.Elements[2], matrix.Elements[3], matrix.OffsetY)
{
}
/// <summary>
/// Creates an instance of this class
/// </summary>
/// <param name="affineTransformVector"></param>
public AffineCoordinateTransformation2D(params double[] affineTransformVector)
{
_a = affineTransformVector[0];
_b = affineTransformVector[1];
_c = affineTransformVector[2];
_d = affineTransformVector[3];
_e = affineTransformVector[4];
_f = affineTransformVector[5];
//var mat = new System.Drawing.Drawing2D.Matrix((float)_a, (float)_b,
// (float)_d, (float)_e, (float)_c, (float)_f);
//var matInv = mat.Clone();
//matInv.Invert();
////
//var m = ToMatrix();
var minv = ToMatrix();
var rep = new alglib.matinv.matinvreport();
var info = 0;
alglib.matinv.rmatrixinverse(ref minv, 3, ref info, rep);
_ainv = minv[0, 0];
_binv = minv[0, 1];
_cinv = minv[2, 0];
_dinv = minv[1, 0];
_einv = minv[1, 1];
_finv = minv[2, 1];
}
public AffineCoordinateTransformation2D(double[,] matrix)
:this(matrix[0,0], matrix[0,1], matrix[2,0],
matrix[1,0], matrix[1,1], matrix[2,1])
{
}
private double[,] ToMatrix()
{
var elem = new double[3,3];
elem[0, 0] = _a;
elem[0, 1] = _b;
//elem[0, 2] = 1d;
elem[1, 0] = _d;
elem[1, 1] = _e;
//elem[1, 2] = 1d;
elem[2, 0] = _c;
elem[2, 1] = _f;
elem[2, 2] = 1d;
return elem;
}
#region Overrides of MathTransform
public override GeoAPI.CoordinateSystems.Transformations.IMathTransform Inverse()
{
//#warning(System.Drawing.Drawing2D.Matrix uses single precision floating point numbers. This involves reduction of precision, not at all accurate!)
//var m = ToMatrix();
//var rep = new alglib.matinv.matinvreport();
//var info = 0;
//alglib.matinv.rmatrixinverse(ref m, 3, ref info, rep);
return new AffineCoordinateTransformation2D(_ainv, _binv, _cinv, _dinv, _einv, _finv);
}
public override double[] Transform(double[] point)
{
/*
Converting from input (X,Y) to output coordinate system (X',Y') is done by:
X' = a*X + b*Y + c,
Y' = d*X + e*Y + f
*/
System.Diagnostics .Debug.Assert(point.Length == 2);
if (IsInverse)
{
return new []
{
_ainv*point[0] + _binv*point[1] + _cinv,
_dinv*point[0] + _einv*point[1] + _finv
};
}
return new []
{
_a*point[0] + _b*point[1] + _c,
_d*point[0] + _e*point[1] + _f
};
}
public override System.Collections.Generic.IList<double[]> TransformList(System.Collections.Generic.IList<double[]> points)
{
System.Collections.Generic.List<System.Double[]> ret = new System.Collections.Generic.List<System.Double[]>(points.Count);
foreach (var d in points)
ret.Add(Transform(d));
return ret;
}
/// <summary>
/// Gets a value indicating whether this transfrom is in inverse mode or not
/// </summary>
public bool IsInverse { get; private set;}
public override void Invert()
{
IsInverse = !IsInverse;
}
public override int DimSource
{
get { throw new System.NotImplementedException(); }
}
public override int DimTarget
{
get { throw new System.NotImplementedException(); }
}
public override string WKT
{
get { throw new System.NotImplementedException(); }
}
public override string XML
{
get { throw new System.NotImplementedException(); }
}
#endregion
#region Implementation of ICoordinateTransformation
public string AreaOfUse
{
get { throw new System.NotImplementedException(); }
}
public string Authority
{
get { throw new System.NotImplementedException(); }
}
public long AuthorityCode
{
get { throw new System.NotImplementedException(); }
}
public GeoAPI.CoordinateSystems.Transformations.IMathTransform MathTransform
{
get { return this; }
}
public string Name
{
get { return "AffineCoordinateTransformation2D"; }
}
public string Remarks
{
get { return ""; }
}
public GeoAPI.CoordinateSystems.ICoordinateSystem SourceCS
{
get { return null; }
}
public GeoAPI.CoordinateSystems.ICoordinateSystem TargetCS
{
get { return null; }
}
public GeoAPI.CoordinateSystems.Transformations.TransformType TransformType
{
get { return GeoAPI.CoordinateSystems.Transformations.TransformType.Transformation; }
}
#endregion
}
[NUnit.Framework.TestFixture]
public class AffineTransformSample
{
private static readonly System.Random RandomNumberGenerator = new System.Random(778564);
static System.Double[] GetRandomOrdinates(System.Int32 size, System.Double min, System.Double max)
{
System.Double[] arr = new System.Double[size];
System.Double width = max - min;
for (System.Int32 i = 0; i < size; i++)
{
System.Double randomValue = RandomNumberGenerator.NextDouble();
arr[i] = min + randomValue * width;
}
return arr;
}
private static SharpMap.Data.FeatureDataTable TransformedFeatureDataTable(
System.Drawing.Drawing2D.Matrix matrix, SharpMap.Data.FeatureDataTable fdt)
{
SharpMap.Data.FeatureDataTable fdtClone = new SharpMap.Data.FeatureDataTable(fdt);
fdtClone.Clear();
foreach (SharpMap.Data.FeatureDataRow row in fdt)
{
SharpMap.Data.FeatureDataRow newRow = fdtClone.NewRow();
for (System.Int32 i = 0; i < fdtClone.Columns.Count; i++)
newRow[i] = row[i];
GeoAPI.Geometries.IPoint smpt = (GeoAPI.Geometries.IPoint)row.Geometry;
System.Drawing.PointF[] pts = new System.Drawing.PointF[]
{ new System.Drawing.PointF((float)smpt.X, (float)smpt.Y) };
matrix.TransformPoints(pts);
newRow.Geometry = new NetTopologySuite.Geometries.Point(pts[0].X, pts[0].Y);
fdtClone.AddRow(newRow);
}
return fdtClone;
}
[NUnit.Framework.Test]
public void TestMatrix()
{
var p = new NetTopologySuite.Geometries.Point(10, 10);
var b = p.AsBinary();
System.Drawing.Drawing2D.Matrix mat = new System.Drawing.Drawing2D.Matrix();
mat.Rotate(30);
mat.Translate(-20, 20);
System.Drawing.PointF[] pts = new System.Drawing.PointF[] { new System.Drawing.PointF(50, 50) };
mat.TransformPoints(pts);
System.Diagnostics.Debug.WriteLine(string.Format("POINT ({0} {1})", pts[0].X, pts[0].Y));
System.Drawing.PointF ptt = pts[0];
System.Drawing.PointF[] ptts = new System.Drawing.PointF[] { new System.Drawing.PointF(ptt.X, ptt.Y) };
System.Drawing.Drawing2D.Matrix inv = mat.Clone();
inv.Invert();
inv.TransformPoints(ptts);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(ptts[0].X - 50f), 0.01);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(ptts[0].Y - 50f), 0.01);
}
[NUnit.Framework.Test]
public void TestMatrix2()
{
System.Drawing.Drawing2D.Matrix mat = new System.Drawing.Drawing2D.Matrix();
mat.Rotate(30);
mat.Translate(-20, 20);
var at = new AffineCoordinateTransformation2D(mat);
var atInv = at.Inverse();
var p0 = new double[] { 50d, 50d };
var pt = at.Transform(p0);
at.Invert();
var p1 = at.Transform(pt);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(p1[0] - p0[0]), 0.01d);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(p1[1] - p0[1]), 0.01d);
var p2 = atInv.Transform(pt);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(p2[0] - p0[0]), 0.01d);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(p2[1] - p0[1]), 0.01d);
System.Drawing.PointF[] pts = new System.Drawing.PointF[] { new System.Drawing.PointF(50, 50) };
mat.TransformPoints(pts);
System.Diagnostics.Debug.WriteLine(string.Format("POINT ({0} {1})", pts[0].X, pts[0].Y));
System.Drawing.PointF ptt = pts[0];
System.Drawing.PointF[] ptts = new System.Drawing.PointF[] { new System.Drawing.PointF(ptt.X, ptt.Y) };
System.Drawing.Drawing2D.Matrix inv = mat.Clone();
inv.Invert();
inv.TransformPoints(ptts);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(ptts[0].X - 50f), 0.01);
NUnit.Framework.Assert.LessOrEqual(System.Math.Abs(ptts[0].Y - 50f), 0.01);
}
[NUnit.Framework.Test]
public void TestAffineTransform2D()
{
//Setup some affine transformation
System.Drawing.Drawing2D.Matrix matrix = new System.Drawing.Drawing2D.Matrix();
matrix.RotateAt(30, new System.Drawing.PointF(0, 0));
matrix.Translate(-20, -20, System.Drawing.Drawing2D.MatrixOrder.Append);
matrix.Shear(0.95f, -0.2f, System.Drawing.Drawing2D.MatrixOrder.Append);
//Create some random sample data
SharpMap.Data.FeatureDataTable fdt1 =
CreatingData.CreatePointFeatureDataTableFromArrays(GetRandomOrdinates(80, -180, 180),
GetRandomOrdinates(80, -90, 90), null);
//Clone random sample data and apply affine transformation on it
SharpMap.Data.FeatureDataTable fdt2 = TransformedFeatureDataTable(matrix, fdt1);
//Get affine transformation with LeastSquaresTransform
SharpMap.Utilities.LeastSquaresTransform lst = new SharpMap.Utilities.LeastSquaresTransform();
//Add at least three corresponding points
lst.AddInputOutputPoint(
((SharpMap.Data.FeatureDataRow)fdt1.Rows[0]).Geometry.Coordinate,
((SharpMap.Data.FeatureDataRow)fdt2.Rows[0]).Geometry.Coordinate);
lst.AddInputOutputPoint(
((SharpMap.Data.FeatureDataRow)fdt1.Rows[39]).Geometry.Coordinate,
((SharpMap.Data.FeatureDataRow)fdt2.Rows[39]).Geometry.Coordinate);
lst.AddInputOutputPoint(
((SharpMap.Data.FeatureDataRow)fdt1.Rows[79]).Geometry.Coordinate,
((SharpMap.Data.FeatureDataRow)fdt2.Rows[79]).Geometry.Coordinate);
/*
//Get affine transformation calculates mean points to improve accuaracy
//Unfortunately the result is not very good, so, since I know better I manually set these
//mean points.
lst.SetMeanPoints(new GeoAPI.Geometries.IPoint(0, 0),
new GeoAPI.Geometries.IPoint(matrix.OffsetX, matrix.OffsetY));
*/
//Create Affine
AffineCoordinateTransformation2D at2 = new AffineCoordinateTransformation2D(lst.GetAffineTransformation());
//Create Map
SharpMap.Map map = new SharpMap.Map(new System.Drawing.Size(720, 360));
//Add not transformed layer
map.Layers.Add(new SharpMap.Layers.VectorLayer("L1",
new SharpMap.Data.Providers.GeometryFeatureProvider(fdt1)));
((SharpMap.Layers.VectorLayer) map.Layers[0]).Style.Symbol =
new System.Drawing.Bitmap(@"..\..\..\DemoWinForm\Resources\flag.png");
//Add transformed layer
map.Layers.Add(new SharpMap.Layers.VectorLayer("L2",
new SharpMap.Data.Providers.GeometryFeatureProvider(fdt2)));
((SharpMap.Layers.VectorLayer) map.Layers[1]).Style.Symbol =
new System.Drawing.Bitmap(@"..\..\..\DemoWinForm\Resources\women.png");
//Render map
map.ZoomToExtents();
//Get map and save to file
var bmp = (System.Drawing.Bitmap)map.GetMap();
bmp.Save("affinetransform1.bmp");
//we want to reverse the previously applied transformation.
((SharpMap.Layers.VectorLayer) map.Layers[1]).CoordinateTransformation = (AffineCoordinateTransformation2D)at2.Inverse();
//Render map
map.ZoomToExtents();
//Get map and save to file
bmp = (System.Drawing.Bitmap)map.GetMap();
bmp.Save("affinetransform2.bmp");
//Hopefully women cover flags ;-).
}
}
#endif
[NUnit.Framework.Test]
public void TestGdalRasterLayer()
{
if (!System.IO.File.Exists("D:\\Daten\\zone49_mga.ecw"))
NUnit.Framework.Assert.Ignore("Adjust file path");
if (!System.IO.File.Exists("D:\\Daten\\zone50_mga.ecw"))
NUnit.Framework.Assert.Ignore("Adjust file path");
var ecw1 = new SharpMap.Layers.GdalRasterLayer("zone49", "D:\\Daten\\zone49_mga.ecw");
var ecw2 = new SharpMap.Layers.GdalRasterLayer("zone50", "D:\\Daten\\zone50_mga.ecw");
var p1 = ecw1.GetProjection();
ecw2.ReprojectToCoordinateSystem(p1);
var m = new SharpMap.Map(new System.Drawing.Size( 1024, 768));
m.Layers.Add(ecw1);
m.Layers.Add(ecw2);
m.ZoomToExtents();
using (var img = m.GetMap())
{
img.Save("ecw.png");
}
}
public static void ReprojectFeatureDataSet(SharpMap.Data.FeatureDataSet fds,
GeoAPI.CoordinateSystems.ICoordinateSystem target)
{
for (var i = 0; i < fds.Tables.Count; i ++)
{
var fdt = fds.Tables[i];
ReprojectFeatureDataTable(fdt, target);
}
}
public static void ReprojectFeatureDataTable(SharpMap.Data.FeatureDataTable fdt,
GeoAPI.CoordinateSystems.ICoordinateSystem target)
{
var source = SharpMap.CoordinateSystems.CoordinateSystemExtensions.GetCoordinateSystem(fdt[0].Geometry);
var ctFactory = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory();
var ct = ctFactory.CreateFromCoordinateSystems(source, target);
var geomFactory = GeoAPI.GeometryServiceProvider.Instance.CreateGeometryFactory((int)target.AuthorityCode);
for (var i = 0; i < fdt.Rows.Count; i++)
{
var fdr = fdt[i];
fdr.Geometry =
GeoAPI.CoordinateSystems.Transformations.GeometryTransform.TransformGeometry(fdr.Geometry,
ct.MathTransform, geomFactory);
}
}
}
}
You can’t perform that action at this time.