言成言成啊 | Kit Chen's Blog

Gis基础知识

发布于2022-05-22 15:13:33，更新于2023-01-30 21:14:16，标签：gis  转载随意，文章会持续修订，请注明来源地址：https://meethigher.top/blog

参考

  1. GIS基础知识 - 坐标系、投影、EPSG:4326、EPSG:3857
  2. GIS基础教程之坐标系 - 知乎
  3. Axis Order — GeoTools 28-SNAPSHOT User Guide
  4. EPSG.io: Coordinate Systems Worldwide
  5. 区域面积-距离/面积计算-示例中心-JS API 示例 | 高德地图API

一、坐标系

1.1 分类

坐标系分为两种

  1. 地理坐标系（Geographic Coordinate System, GCS）
  2. 投影坐标系（Projected Coordinate System, PCS）

这边知乎老哥的文章描述很容易理解，我直接cv过来。

地理坐标系，就是一个球，说成球面坐标系反而好理解一点。就比如一个完整的橘子。

投影坐标系，就是将球展开后的一个面。就比如橘子剥皮展开。

地理坐标系和投影坐标系的最大区分，就是单位不同。

地理坐标系以度为单位。

投影坐标系以米为单位。

1.2 常用坐标系

对于地图开发人员来说，常见的测量系统就是WGS84测量系统(1984年全球测量系统)。

基于WGS84的常用坐标系为两种。

4326是地理坐标系，以度为单位。广泛用于GPS定位等。

3857是投影坐标系，以米为单位。用于平面地图，如GoogleMap、OpenStreetMap、BingMap等。

如果想要找取专属于自己地区、更精确的坐标系。可以去epsg.io搜索。

二、geotools

作为用java做饭碗的人来说，想要研究Gis，肯定少不了跟geotools打交道。

2.1 maven依赖

pom.xml

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
<repositories>
    <repository>
        <id>osgeo</id>
        <name>OSGeo Release Repository</name>
        <url>https://repo.osgeo.org/repository/release/</url>
        <snapshots>
            <enabled>false</enabled>
        </snapshots>
        <releases>
            <enabled>true</enabled>
        </releases>
    </repository>
</repositories>

<dependencies>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-grid</artifactId>
        <version>${geotools.version}</version>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-opengis</artifactId>
        <version>${geotools.version}</version>
    </dependency>


    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-epsg-hsql</artifactId>
        <version>${geotools.version}</version>
        <exclusions>
            <exclusion>
                <groupId>javax</groupId>
                <artifactId>javaee-api</artifactId>
            </exclusion>
        </exclusions>
    </dependency>
    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-main</artifactId>
        <version>${geotools.version}</version>
        <exclusions>
            <exclusion>
                <groupId>javax</groupId>
                <artifactId>javaee-api</artifactId>
            </exclusion>
        </exclusions>
    </dependency>

    <dependency>
        <groupId>org.geotools</groupId>
        <artifactId>gt-geojson</artifactId>
        <version>${geotools.version}</version>
        <exclusions>
            <exclusion>
                <groupId>javax</groupId>
                <artifactId>javaee-api</artifactId>
            </exclusion>
        </exclusions>
    </dependency>

    <dependency>
        <groupId>junit</groupId>
        <artifactId>junit</artifactId>
        <version>4.12</version>
        <scope>test</scope>
    </dependency>
</dependencies>

2.2 简单使用

下面简单演示一些常见用法。我这边就全以wkt为例了。wkt的绘制可以参考wkt在线绘制展示_EPSG4326。

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
/**
 * GeoTools简单使用
 *
 * @author chenchuancheng github.com/meethigher
 * @since 2022/5/22 20:30
 */
public class GeotoolsUsing {

    @Test
    public void testCreateGeo() throws Exception {
        //这是4326坐标系的wkt
        String wkt = "LINESTRING(110.26968650499015 35.774090269456565,118.2156404214427 35.973104738241815)";
        WKTReader reader = new WKTReader();
        Geometry geometry = reader.read(wkt);
        System.out.println(geometry.buffer(1));
        Assert.assertNotNull(geometry);
    }

    @Test
    public void testGeoApi() throws Exception {
        //area1是包含area2的
        String area1 = "POLYGON((113.36572265625 38.80204526520603,107.28369301557538 33.994384090476174,112.38134926557541 28.960087983961387,119.2368153333664 31.690781121452815,113.36572265625 38.80204526520603))";
        String area2 = "LINESTRING(110.35095139165554 33.90070434182671,116.21651043156132 33.564288629604945)";
        WKTReader reader = new WKTReader();
        Geometry geometry1 = reader.read(area1);
        Geometry geometry2 = reader.read(area2);

        System.out.println(geometry1.contains(geometry2));//true

        //获取边界，即获取最北、最南、最东、最西，组成一个矩形
        Geometry envelope = geometry1.getEnvelope();
        System.out.println(envelope.toString());//输出wkt

        //获取边界的点
        Envelope internal = geometry1.getEnvelopeInternal();
        double minX = internal.getMinX();
        double minY = internal.getMinY();

        //返回其中一个顶点
        Coordinate coordinate = geometry1.getCoordinate();
        System.out.println(coordinate);

        //拓宽区域，由于使用的是4326坐标系，1在此处表示1度。这个单位也是跟着坐标系走的
        //以面的边向外拓宽1，顶点处就是以1为半径的圆弧。如果是线，就是四周拓宽
        Geometry buffer = geometry1.buffer(1);
        System.out.println(buffer);

        //判断交叉，我画的是个不相交，但包含的。
        //这个方法有两个含义。a.相交 b. 包含 满足其一即可
        boolean intersects = geometry1.intersects(geometry2);
        System.out.println(intersects);//true

        //面求面积，线求长度。单位同样是跟着坐标系走
        System.out.println(geometry1.getArea());
        System.out.println(geometry2.getLength());
    }


    @Test
    public void testTransformCrs() throws Exception {
        //我使用的wkt都是经度在前，一般正常使用也是经度在前。将4326转换为3857计算面积(㎡)
        CRSAuthorityFactory factory = CRS.getAuthorityFactory(true);
        CoordinateReferenceSystem source = factory.createCoordinateReferenceSystem("EPSG:4326");
        CoordinateReferenceSystem target = factory.createCoordinateReferenceSystem("EPSG:3857");
        MathTransform transform = CRS.findMathTransform(source, target);

        //高德地图面积https://developer.amap.com/demo/javascript-api/example/calcutation/ring-area
        // 该区域高德面积125178908819平方米，坐标系不同，有误差很正常
        String wkt = "POLYGON((110.4141583943802 36.03154311769477,110.54418331527785 33.09352034423874,114.95057594167424 33.32318635566811,114.73386810758919 35.99648515314789,110.4141583943802 36.03154311769477))";

        WKTReader reader = new WKTReader();
        Geometry geometry = reader.read(wkt);
        System.out.println(geometry.getArea());//度为单位，12.258110818164612
        Geometry geometry1 = JTS.transform(geometry, transform);
        //用bigdecimal只是单纯不想看到科学计数法。
        System.out.println(BigDecimal.valueOf(geometry1.getArea()));//米为单位，184591949648.2439
    }
}

三、常用坐标系转换

3.1 java工具类

抄袭自JAVA实现WGS84、百度坐标系、高德坐标系转化工具类_liuccn的博客-CSDN博客_java百度坐标转高德坐标

进行修改后如下

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
/**
 * 仅适用于中国的点位坐标系转换
 * 各地图API坐标系统比较与转换;
 * <p>
 * WGS84坐标系(EPSG4326)：即地球坐标系，国际上通用的坐标系。设备一般包含GPS芯片或者北斗芯片获取的经纬度为WGS84地理坐标系,
 * 谷歌地图采用的是WGS84地理坐标系（中国范围除外）;
 * <p>
 * GCJ02坐标系：即火星坐标系，是由中国国家测绘局制订的地理信息系统的坐标系统。由WGS84坐标系经加密后的坐标系。
 * 高德地图、腾讯地图、谷歌中国地图和搜搜中国地图采用的是GCJ02地理坐标系;
 * <p>
 * BD09坐标系：即百度坐标系，GCJ02坐标系经加密后的坐标系;
 * <p>
 * 搜狗坐标系、图吧坐标系等，估计也是在GCJ02基础上加密而成的。
 *
 * @author https://blog.csdn.net/lc_2014c/article/details/125878730
 * @copyer chenchuancheng
 * @since 2023/1/6 16:13
 */
public strictfp class PointConvert {

    //private final static double pi = 3.1415926535897932384626;
    private final static double pi = Math.PI;
    private final static double a = 6378245.0;
    private final static double ee = 0.00669342162296594323;

    /**
     * 是否超出中国范围以外
     *
     * @param lat 纬度
     * @param lon 经度
     * @return true表示中国外
     */
    private static boolean outOfChina(double lat, double lon) {
        if (lon < 72.004 || lon > 137.8347)
            return true;
        return lat < 0.8293 || lat > 55.8271;
    }


    /**
     * 转换纬度
     *
     * @param x
     * @param y
     * @return
     */
    private static double transformLat(double x, double y) {
        double ret = -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * x * y
                + 0.2 * Math.sqrt(Math.abs(x));
        ret += (20.0 * Math.sin(6.0 * x * pi) + 20.0 * Math.sin(2.0 * x * pi)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(y * pi) + 40.0 * Math.sin(y / 3.0 * pi)) * 2.0 / 3.0;
        ret += (160.0 * Math.sin(y / 12.0 * pi) + 320 * Math.sin(y * pi / 30.0)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * 转换经度
     *
     * @param x
     * @param y
     * @return
     */
    private static double transformLon(double x, double y) {
        double ret = 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * x * y + 0.1
                * Math.sqrt(Math.abs(x));
        ret += (20.0 * Math.sin(6.0 * x * pi) + 20.0 * Math.sin(2.0 * x * pi)) * 2.0 / 3.0;
        ret += (20.0 * Math.sin(x * pi) + 40.0 * Math.sin(x / 3.0 * pi)) * 2.0 / 3.0;
        ret += (150.0 * Math.sin(x / 12.0 * pi) + 300.0 * Math.sin(x / 30.0
                * pi)) * 2.0 / 3.0;
        return ret;
    }

    /**
     * 转换为84坐标系
     *
     * @param lat
     * @param lon
     * @return
     */
    private static Point transform(double lat, double lon) {
        if (outOfChina(lat, lon)) {
            return new Point(lat, lon);
        }
        double dLat = transformLat(lon - 105.0, lat - 35.0);
        double dLon = transformLon(lon - 105.0, lat - 35.0);
        double radLat = lat / 180.0 * pi;
        double magic = Math.sin(radLat);
        magic = 1 - ee * magic * magic;
        double sqrtMagic = Math.sqrt(magic);
        dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
        dLon = (dLon * 180.0) / (a / sqrtMagic * Math.cos(radLat) * pi);
        double mgLat = lat + dLat;
        double mgLon = lon + dLon;

        return new Point(mgLat, mgLon);
    }


    /**
     * 84 to 火星坐标系 (GCJ-02)
     * <p>
     * World Geodetic System ==> Mars Geodetic System
     *
     * @param lat
     * @param lon
     * @return
     */
    public static Point wgs84_to_gcj02(double lat, double lon) {
        if (outOfChina(lat, lon)) {
            return null;
        }
        double dLat = transformLat(lon - 105.0, lat - 35.0);
        double dLon = transformLon(lon - 105.0, lat - 35.0);
        double radLat = lat / 180.0 * pi;
        double magic = Math.sin(radLat);
        magic = 1 - ee * magic * magic;
        double sqrtMagic = Math.sqrt(magic);
        dLat = (dLat * 180.0) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi);
        dLon = (dLon * 180.0) / (a / sqrtMagic * Math.cos(radLat) * pi);
        double mgLat = lat + dLat;
        double mgLon = lon + dLon;

        return new Point(mgLat, mgLon);
    }

    /**
     * 火星坐标系 (GCJ-02) to 84
     *
     * @param lon
     * @param lat
     * @return
     */
    public static Point gcj02_to_wgs84(double lat, double lon) {
        Point gps = transform(lat, lon);
        double longitude = lon * 2 - gps.getLon();
        double latitude = lat * 2 - gps.getLat();
        return new Point(latitude, longitude);
    }


    /**
     * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法
     * <p>
     * 将 GCJ-02 坐标转换成 BD-09 坐标
     *
     * @param gg_lat
     * @param gg_lon
     */
    public static Point gcj02_to_bd09(double gg_lat, double gg_lon) {
        double x = gg_lon, y = gg_lat;

        double z = Math.sqrt(x * x + y * y) + 0.00002 * Math.sin(y * pi);

        double theta = Math.atan2(y, x) + 0.000003 * Math.cos(x * pi);

        double bd_lon = z * Math.cos(theta) + 0.0065;

        double bd_lat = z * Math.sin(theta) + 0.006;

        return new Point(bd_lat, bd_lon);
    }

    /**
     * 火星坐标系 (GCJ-02) 与百度坐标系 (BD-09) 的转换算法
     * <p>
     * 将 BD-09 坐标转换成GCJ-02 坐标
     *
     * @param bd_lat
     * @param bd_lon
     * @return
     */
    public static Point bd09_to_gcj02(double bd_lat, double bd_lon) {
        double x = bd_lon - 0.0065, y = bd_lat - 0.006;
        double z = Math.sqrt(x * x + y * y) - 0.00002 * Math.sin(y * pi);

        double theta = Math.atan2(y, x) - 0.000003 * Math.cos(x * pi);

        double gg_lon = z * Math.cos(theta);

        double gg_lat = z * Math.sin(theta);

        return new Point(gg_lat, gg_lon);
    }


    /**
     * WGS84坐标系(EPSG:4326)转换为墨卡托投影坐标系(EPSG:3857)
     */
    public static Point wgs84_to_mercator(double lon, double lat) {
        //wgs84下的地球半径
        double earthRadius = 6378137.0;
        double arc = lat * pi / 180;
        return new Point(lon * pi / 180 * earthRadius,
                earthRadius / 2 * Math.log((1.0 + Math.sin(arc)) / (1.0 - Math.sin(arc))));
    }

    /**
     * 墨卡托投影坐标系(EPSG:3857)转换为wgs84坐标系(EPSG:4326)
     */
    public static Point mercator_to_wgs84(double lon, double lat) {
        //wgs84下的地球半径
        double earthRadius = 6378137.0;
        //周长
        double round = 2 * pi * earthRadius / 2.0;
        return new Point(lat / round * 180,
                180 / pi * (2 * Math.atan(Math.exp((lon / round * 180) * Math.PI / 180)) - Math.PI / 2));
    }

    /**
     * 点位信息
     *
     * @author chenchuancheng
     * @since 2023/1/6 16:14
     */
    public final static class Point {

        /**
         * 经度
         */
        private double lon;

        /**
         * 纬度
         */
        private double lat;

        public Point(double lat, double lon) {
            this.lon = lon;
            this.lat = lat;
        }

        public double getLon() {
            return lon;
        }

        public void setLon(double lon) {
            this.lon = lon;
        }

        public double getLat() {
            return lat;
        }

        public void setLat(double lat) {
            this.lat = lat;
        }
    }
}

以火星坐标系GCJ02转换WGS84为例，准确性如下。

3.2 postgis方法

抄袭自geocompass/pg-coordtransform: 基于PostgreSQL+PostGIS的火星坐标系、百度坐标系、WGS84坐标系、CGCS2000坐标系的转换函数

开启postgis扩展

1
create extension postgis

创建方法

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
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
1001
1002
1003
1004
1005
1006
1007
1008
1009
1010
1011
1012
1013
1014
1015
1016
1017
1018
1019
1020
1021
1022
1023
1024
1025
1026
1027
1028
1029
1030
1031
1032
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043
1044
1045
1046
1047
1048
1049
1050
1051
1052
1053
1054
1055
1056
1057
1058
1059
1060
1061
1062
1063
1064
1065
1066
1067
1068
1069
1070
1071
1072
1073
1074
1075
1076
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090
1091
1092
1093
1094
1095
1096
1097
1098
1099
1100
1101
1102
1103
1104
1105
1106
1107
1108
1109
1110
1111
1112
1113
1114
1115
1116
1117
1118
1119
1120
1121
1122
1123
1124
1125
CREATE OR REPLACE FUNCTION "public"."geoc_bd09towgs84_multipolygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    single_polygon  geometry;
    single_polygon_trans  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_MultiPolygon' THEN
        RETURN null;
    END IF;
    FOR single_polygon IN SELECT (ST_Dump($1)).geom LOOP
        single_polygon_trans := geoc_bd09towgs84_polygon(single_polygon); 
        target_parts := array_append(target_parts,single_polygon_trans);
    END LOOP;
    
    SELECT st_multi(ST_Union(target_parts)) INTO final_geom;
    RETURN final_geom;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09towgs84_point"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
  x       numeric; 
	y       numeric; 
	gcj_point  geometry;
	wgs_point  geometry;

BEGIN
    if st_geometrytype(geom) != 'ST_Point' then
      return null;
    end if;
    x := st_x(geom);
    y := st_y(geom);
    if (geoc_is_in_china_bbox(x, y) = false) then
      return geom;
    end if;
    gcj_point = geoc_bd09togcj02_point(geom);
	  wgs_point = geoc_gcj02towgs84_point(gcj_point);
	  return wgs_point;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09towgs84_polygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    source_npoints  integer;
    single_line     geometry;
    single_line_trans geometry;
    single_polygon  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_Polygon' THEN
        RETURN null;
    END IF;

    FOR single_polygon IN SELECT ST_ExteriorRing ((st_dumprings($1)).geom) as geom LOOP				
        source_npoints := ST_NPoints(single_polygon); 
        single_line  := ST_RemovePoint(single_polygon, source_npoints - 1);  
        single_line_trans := geoc_bd09towgs84_line(single_line);  
        target_parts := array_append(target_parts, ST_AddPoint(single_line_trans, ST_PointN(single_line_trans, 1)));  
    END LOOP;
    SELECT ST_MakePolygon(target_parts[1], target_parts[2:array_upper(target_parts, 1)]) INTO final_geom; 
    RETURN final_geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_cgcs2000tobd09"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
BEGIN
IF st_srid(geom) != '4490' THEN
        RETURN null;
end if;
return geoc_wgs84tobd09(st_transform(st_setsrid(geom,4490),4326));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_cgcs2000togcj02"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
BEGIN
IF st_srid(geom) != '4490' THEN
        RETURN null;
end if;
return geoc_wgs84togcj02(st_transform(st_setsrid(geom,4490),4326));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_delta"("lon" numeric, "lat" numeric)
  RETURNS "pg_catalog"."jsonb" AS $BODY$
DECLARE
    ret             varchar;
    dLon            numeric;
    dlat            numeric;
    radLat          numeric;
    magic           numeric;
    sqrtMagic       numeric;
    ee              numeric;
    a               numeric;
BEGIN
    ee := 0.006693421622965823;
    a  := 6378245;
    dLon := geoc_transform_lon(lon - 105, lat - 35);
    dLat := geoc_transform_lat(lon - 105, lat - 35);
    radLat := lat / 180 * pi();
    magic = sin(radLat);

    magic = 1 - ee * magic * magic;

    sqrtMagic := sqrt(magic);
    dLon = (dLon * 180) / (a / sqrtMagic * cos(radLat) * pi());
    dLat = (dLat * 180) / ((a * (1 - ee)) / (magic * sqrtMagic) * pi());
    ret :='['||dLon||','||dLat||']';
	return ret::jsonb;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tobd09"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
BEGIN
IF st_srid(geom) != '4490' and  st_srid(geom) != '4326'THEN
        RETURN null;
end if;
case ST_GeometryType(geom)
    when 'ST_LineString' then 
			return geoc_gcj02tobd09_line(geom);
	when 'ST_MultiLineString' then 
		return geoc_gcj02tobd09_multiline(geom);
	when 'ST_Point' then 
		return geoc_gcj02tobd09_point(geom);
	when 'ST_MultiPoint' then 
			return geoc_gcj02tobd09_multipoint(geom);
	when 'ST_Polygon' then 
			return geoc_gcj02tobd09_polygon(geom);
	when 'ST_MultiPolygon' then
		return geoc_gcj02tobd09_multipolygon(geom); 
	ELSE
    	RETURN null;
END CASE;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tobd09_line"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		p_p     geometry;
		p_t     geometry;
		z_t     geometry;
		i       int;
BEGIN
    i:=1;
	while i <= st_npoints(geom) LOOP
		p_p := st_pointn(geom,i);
		p_t := geoc_gcj02tobd09_point(p_p);
		geom:=st_setpoint(geom,i-1,p_t);
		i:=i+1;
	end LOOP;
	return geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tobd09_multiline"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom  $Q$ using geom LOOP
	  	transform_i :=geoc_gcj02tobd09_line(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tobd09_multipoint"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom  $Q$ using geom LOOP
	  	transform_i :=geoc_gcj02tobd09_point(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tobd09_multipolygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    single_polygon  geometry;
    single_polygon_trans  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_MultiPolygon' THEN
        RETURN null;
    END IF;
        FOR single_polygon IN SELECT (ST_Dump($1)).geom LOOP
                single_polygon_trans := geoc_gcj02tobd09_polygon(single_polygon); 
                target_parts := array_append(target_parts,single_polygon_trans);
        END LOOP;
				
        SELECT st_multi(ST_Union(target_parts)) INTO final_geom;
        RETURN final_geom;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tobd09_point"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
z         double precision; 
theta     double precision; 
x_pi      double precision:=3.14159265358979324 * 3000.0 / 180.0; 
lon       numeric; 
lat       numeric; 
bd_point  geometry;

BEGIN
    if st_geometrytype(geom) != 'ST_Point' then
      return null;
    end if;
    lon := st_x(geom);
    lat := st_y(geom);
    if geoc_is_in_china_bbox(lon, lat) = false THEN 
      return geom;
    end if;
    z:= sqrt(power(lon,2) + power(lat,2)) + 0.00002 * sin(lat * x_pi); 
    theta:= atan2(lat, lon) + 0.000003 * cos(lon * x_pi); 
	bd_point:=ST_SetSRID(ST_MakePoint(z * cos(theta) + 0.0065,z * sin(theta) + 0.006),4326); 
	  return bd_point; 

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tobd09_polygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    source_npoints  integer;
    single_line     geometry;
    single_line_trans geometry;
    single_polygon  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_Polygon' THEN
        RETURN null;
    END IF;

    FOR single_polygon IN SELECT ST_ExteriorRing ((st_dumprings($1)).geom) as geom LOOP				
        source_npoints := ST_NPoints(single_polygon); 
        single_line  := ST_RemovePoint(single_polygon, source_npoints - 1);  
        single_line_trans := geoc_gcj02tobd09_line(single_line);  
        target_parts := array_append(target_parts, ST_AddPoint(single_line_trans, ST_PointN(single_line_trans, 1)));  
    END LOOP;
    SELECT ST_MakePolygon(target_parts[1], target_parts[2:array_upper(target_parts, 1)]) INTO final_geom;
    RETURN final_geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02tocgcs2000"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
BEGIN
IF st_srid(geom) != '4490' THEN
        RETURN null;
end if;
return st_transform(st_setsrid(geoc_gcj02towgs84(geom),4326),4490);
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02towgs84"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
BEGIN
IF st_srid(geom) != '4490' and  st_srid(geom) != '4326'THEN
        RETURN null;
end if;
case ST_GeometryType(geom)
    when 'ST_LineString' then 
			return geoc_gcj02towgs84_line(geom);
	when 'ST_MultiLineString' then 
		return geoc_gcj02towgs84_multiline(geom);
	when 'ST_Point' then 
		return geoc_gcj02towgs84_point(geom);
	when 'ST_MultiPoint' then 
			return geoc_gcj02towgs84_multipoint(geom);
	when 'ST_Polygon' then 
			return geoc_gcj02towgs84_polygon(geom);
	when 'ST_MultiPolygon' then
		return geoc_gcj02towgs84_multipolygon(geom); 
	ELSE
    	RETURN null;
END CASE;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02towgs84_line"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		p_p     geometry;
		p_t     geometry;
		z_t     geometry;
		i       int;
BEGIN
    i:=1;
	while i <= st_npoints(geom) LOOP
		p_p := st_pointn(geom,i);
		p_t := geoc_gcj02towgs84_point(p_p);
		geom:=st_setpoint(geom,i-1,p_t);
		i:=i+1;
	end LOOP;
	return geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02towgs84_multiline"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom  $Q$ using geom LOOP
	  	transform_i :=geoc_gcj02towgs84_line(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02towgs84_multipoint"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom  $Q$ using geom LOOP
	  	transform_i :=geoc_gcj02towgs84_point(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02towgs84_multipolygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    single_polygon  geometry;
    single_polygon_trans  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_MultiPolygon' THEN
        RETURN null;
    END IF;
        FOR single_polygon IN SELECT (ST_Dump($1)).geom LOOP
                single_polygon_trans := geoc_gcj02towgs84_polygon(single_polygon); 
                target_parts := array_append(target_parts,single_polygon_trans);
        END LOOP;
				
        SELECT st_multi(ST_Union(target_parts)) INTO final_geom;
        RETURN final_geom;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02towgs84_point"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    tempPoint  numeric[];
    wgsLon     numeric;
    wgsLat     numeric;
    lon        numeric;
    lat        numeric;
BEGIN
    if st_geometrytype(geom) != 'ST_Point' then
      return null;
    end if;
    lon := st_x(geom);
    lat := st_y(geom);
    if geoc_is_in_china_bbox(lon, lat) = false THEN 
      return geom;
    end if;

    tempPoint := geoc_wgs84togcj02(ARRAY[lon, lat]);
    wgsLon := lon*2 - tempPoint[1];
    wgsLat := lat*2 - tempPoint[2];
    return st_makepoint(wgsLon,wgsLat);
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_gcj02towgs84_polygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    source_npoints  integer;
    single_line     geometry;
    single_line_trans geometry;
    single_polygon  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_Polygon' THEN
        RETURN null;
    END IF;

    FOR single_polygon IN SELECT ST_ExteriorRing ((st_dumprings($1)).geom) as geom LOOP				
        source_npoints := ST_NPoints(single_polygon); 
        single_line  := ST_RemovePoint(single_polygon, source_npoints - 1);  
        single_line_trans := geoc_gcj02towgs84_line(single_line);  
        target_parts := array_append(target_parts, ST_AddPoint(single_line_trans, ST_PointN(single_line_trans, 1)));  
    END LOOP;
    SELECT ST_MakePolygon(target_parts[1], target_parts[2:array_upper(target_parts, 1)]) INTO final_geom; 
    RETURN final_geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_is_in_china_bbox"("lon" numeric, "lat" numeric)
  RETURNS "pg_catalog"."bool" AS $BODY$
DECLARE
BEGIN
    
    return lon >= 72.004 and lon <= 137.8347 and lat >= 0.8293 and lat <= 55.8271;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_transform_lat"("x" numeric, "y" numeric)
  RETURNS "pg_catalog"."numeric" AS $BODY$
DECLARE
ret   numeric;
BEGIN
    ret := -100 + 2 * x + 3 * y + 0.2 * y * y + 0.1 * x * y + 0.2 * sqrt(abs(x));
    ret := ret + (20 * sin(6 * x * PI()) + 20 * sin(2 * x * PI())) * 2 / 3;
    ret := ret +(20 * sin(y * PI()) + 40 * sin(y / 3 * PI())) * 2 / 3;
    ret := ret +(160 * sin(y / 12 * PI()) + 320 * sin(y * PI() / 30)) * 2 / 3;
    return ret;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_transform_lon"("x" numeric, "y" numeric)
  RETURNS "pg_catalog"."numeric" AS $BODY$
DECLARE
    ret   numeric;
BEGIN
    ret := 300 + x + 2 * y + 0.1 * x * x + 0.1 * x * y + 0.1 * sqrt(abs(x));
    ret :=ret + (20 * sin(6 * x * pi()) + 20 * sin(2 * x * pi())) * 2 / 3;
    ret :=ret + (20 * sin(x * pi()) + 40 * sin(x / 3 * pi())) * 2 / 3;
    ret :=ret + (150 * sin(x / 12 * pi()) + 300 * sin(x / 30 * pi())) * 2 / 3;
    return ret;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84tobd09"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
BEGIN
IF st_srid(geom) != '4490' and  st_srid(geom) != '4326'THEN
        RETURN null;
end if;
case ST_GeometryType(geom)
  when 'ST_LineString' then 
		return geoc_wgs84tobd09_line(geom);
	when 'ST_MultiLineString' then 
		return geoc_wgs84tobd09_multiline(geom);
	when 'ST_Point' then 
		return geoc_wgs84tobd09_point(geom);
	when 'ST_MultiPoint' then 
		return geoc_wgs84tobd09_multipoint(geom);
	when 'ST_Polygon' then 
		return geoc_wgs84tobd09_polygon(geom);
	when 'ST_MultiPolygon' then
		return geoc_wgs84tobd09_multipolygon(geom); 
	ELSE
    	RETURN null;
END CASE;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84tobd09_line"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		p_p     geometry;
		p_t     geometry;
		z_t     geometry;
		i       int;
BEGIN
    i:=1;
	while i <= st_npoints(geom) LOOP
		p_p := st_pointn(geom,i);
		p_t := geoc_wgs84tobd09_point(p_p);
		geom:=st_setpoint(geom,i-1,p_t);
		i:=i+1;
	end LOOP;
	return geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84tobd09_multiline"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom  $Q$ using geom LOOP
	  	transform_i :=geoc_wgs84tobd09_line(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84tobd09_multipoint"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom  $Q$ using geom LOOP
	  	transform_i :=geoc_wgs84tobd09_point(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84tobd09_multipolygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    single_polygon  geometry;
    single_polygon_trans  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_MultiPolygon' THEN
        RETURN null;
    END IF;
        FOR single_polygon IN SELECT (ST_Dump($1)).geom LOOP
                single_polygon_trans := geoc_wgs84tobd09_polygon(single_polygon); 
                target_parts := array_append(target_parts,single_polygon_trans);
        END LOOP;
				
        SELECT st_multi(ST_Union(target_parts)) INTO final_geom;
        RETURN final_geom;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84tobd09_point"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
lon       numeric; 
lat       numeric; 
bd_point  geometry;
gcj_point  geometry;

BEGIN
    if st_geometrytype(geom) != 'ST_Point' then
      return null;
    end if;
    lon := st_x(geom);
    lat := st_y(geom);
    if geoc_is_in_china_bbox(lon, lat) = false THEN 
      return geom;
    end if;
    gcj_point = geoc_wgs84togcj02_point(geom);
		bd_point = geoc_gcj02tobd09_point(gcj_point);
	  return bd_point; 

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84tobd09_polygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    source_npoints  integer;
    single_line     geometry;
    single_line_trans geometry;
    single_polygon  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_Polygon' THEN
        RETURN null;
    END IF;

    FOR single_polygon IN SELECT ST_ExteriorRing ((st_dumprings($1)).geom) as geom LOOP				
        source_npoints := ST_NPoints(single_polygon); 
        single_line  := ST_RemovePoint(single_polygon, source_npoints - 1);  
        single_line_trans := geoc_wgs84tobd09_line(single_line);  
        target_parts := array_append(target_parts, ST_AddPoint(single_line_trans, ST_PointN(single_line_trans, 1)));  
    END LOOP;
    SELECT ST_MakePolygon(target_parts[1], target_parts[2:array_upper(target_parts, 1)]) INTO final_geom;  
    RETURN final_geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
IF st_srid(geom) != '4490' and  st_srid(geom) != '4326'THEN
        RETURN null;
end if;
 	CASE ST_GeometryType(geom)
    	when 'ST_LineString' then 
			return geoc_wgs84togcj02_line(geom);
		when 'ST_MultiLineString' then 
			return geoc_wgs84togcj02_multiline(geom);
		when 'ST_Point' then 
			return geoc_wgs84togcj02_point(geom);
		when 'ST_MultiPoint' then 
			return geoc_wgs84togcj02_multipoint(geom);
		when 'ST_Polygon' then 
			return geoc_wgs84togcj02_polygon(geom);
		when 'ST_MultiPolygon' then
			return geoc_wgs84togcj02_multipolygon(geom);
		ELSE
     		RETURN null;
	END CASE;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02_line"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		p_p     geometry;
		p_t     geometry;
		z_t     geometry;
		i       int;
BEGIN
    i:=1;
	while i <= st_npoints(geom) LOOP
		p_p := st_pointn(geom,i);
    	p_t := geoc_wgs84togcj02_point(p_p);
		geom:=st_setpoint(geom,i-1,p_t);
		i:=i+1;
	end LOOP;
	return geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02_multiline"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom $Q$ using geom LOOP
	  	transform_i :=geoc_wgs84togcj02_line(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02_multipoint"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom $Q$ using geom LOOP
	  	transform_i :=geoc_wgs84togcj02_point(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02_multipolygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    single_polygon  geometry;
    single_polygon_trans  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_MultiPolygon' THEN
        RETURN null;
    END IF;
    FOR single_polygon IN SELECT (ST_Dump($1)).geom LOOP
        single_polygon_trans := geoc_wgs84togcj02_polygon(single_polygon); 
        target_parts := array_append(target_parts,single_polygon_trans);
    END LOOP;
    
    SELECT st_multi(ST_Union(target_parts)) INTO final_geom;
    RETURN final_geom;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02"("coord" _numeric)
  RETURNS "pg_catalog"."_numeric" AS $BODY$
DECLARE
    ret             numeric[];
    dLon            numeric;
    dlat            numeric;
		lon             numeric;
		lat             numeric;
		d								jsonb;
-- 		coord           ARRAY;
BEGIN
    lon := coord[1];
    lat := coord[2];
    if (geoc_is_in_china_bbox(lon, lat) = false) then
        return coord;
    end if;
    d := geoc_delta(lon, lat);
    dlon := d->0;
    dlat := d->1;
    ret := ARRAY[lon + dlon , lat + dlat];
    return ret;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02_point"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    lon     numeric;
    lat     numeric;
    d       jsonb;
    dlon    numeric;
    dlat    numeric;
BEGIN
    if st_geometrytype(geom) != 'ST_Point' then
        return null;
    end if;
    lon := st_x(geom);
    lat := st_y(geom);
    if (geoc_is_in_china_bbox(lon, lat) = false) then
        return geom;
    end if;
    d := geoc_delta(lon, lat);
    dlon := d->0;
    dlat := d->1;
    return st_makepoint(lon + dlon,lat + dlat);
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_wgs84togcj02_polygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    source_npoints  integer;
    single_line     geometry;
    single_line_trans geometry;
    single_polygon  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_Polygon' THEN
        RETURN null;
    END IF;

    FOR single_polygon IN SELECT ST_ExteriorRing ((st_dumprings($1)).geom) as geom LOOP				
        source_npoints := ST_NPoints(single_polygon); 
        single_line  := ST_RemovePoint(single_polygon, source_npoints - 1);  
        single_line_trans := geoc_wgs84togcj02_line(single_line);  
        target_parts := array_append(target_parts, ST_AddPoint(single_line_trans, ST_PointN(single_line_trans, 1)));  
    END LOOP;
    SELECT ST_MakePolygon(target_parts[1], target_parts[2:array_upper(target_parts, 1)]) INTO final_geom; 
    RETURN final_geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09tocgcs2000"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
BEGIN
IF st_srid(geom) != '4490' THEN
        RETURN null;
end if;
return st_transform(st_setsrid(geoc_bd09towgs84(geom),4326),4490);
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09togcj02"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
IF st_srid(geom) != '4490' and  st_srid(geom) != '4326'THEN
        RETURN null;
	end if;
 	CASE ST_GeometryType(geom)
    	when 'ST_LineString' then 
			return geoc_bd09togcj02_line(geom);
		when 'ST_MultiLineString' then 
			return geoc_bd09togcj02_multiline(geom);
		when 'ST_Point' then 
			return geoc_bd09togcj02_point(geom);
		when 'ST_MultiPoint' then 
			return geoc_bd09togcj02_multipoint(geom);
		when 'ST_Polygon' then 
			return geoc_bd09togcj02_polygon(geom);
		when 'ST_MultiPolygon' then
			return geoc_bd09togcj02_multipolygon(geom);
		ELSE
     		RETURN null;
	END CASE;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09togcj02_line"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		p_p     geometry;
		p_t     geometry;
		z_t     geometry;
		i       int;
BEGIN
    i:=1;
	while i <= st_npoints(geom) LOOP
		p_p := st_pointn(geom,i);
    	p_t := geoc_bd09togcj02_point(p_p);
		geom:=st_setpoint(geom,i-1,p_t);
		i:=i+1;
	end LOOP;
	return geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09togcj02_multiline"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom $Q$ using geom LOOP
	  	transform_i :=geoc_bd09togcj02_line(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09togcj02_multipoint"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom $Q$ using geom LOOP
	  	transform_i :=geoc_bd09togcj02_point(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09togcj02_multipolygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    single_polygon  geometry;
    single_polygon_trans  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_MultiPolygon' THEN
        RETURN null;
    END IF;
    FOR single_polygon IN SELECT (ST_Dump($1)).geom LOOP
        single_polygon_trans := geoc_bd09togcj02_polygon(single_polygon); 
        target_parts := array_append(target_parts,single_polygon_trans);
    END LOOP;
    
    SELECT st_multi(ST_Union(target_parts)) INTO final_geom;
    
    RETURN final_geom;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09togcj02_point"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
  x       numeric; 
	y       numeric; 
 	z       double precision; 
 	theta   double precision; 
 	x_pi    double precision:=3.14159265358979324 * 3000.0 / 180.0; 
	gcj_point  geometry;

BEGIN
    if st_geometrytype(geom) != 'ST_Point' then
        return null;
    end if;
    x := st_x(geom);
    y := st_y(geom);
    if (geoc_is_in_china_bbox(x, y) = false) then
        return geom;
    end if;
    x:= ST_X(geom) - 0.0065; 
    y:= ST_Y(geom) - 0.006; 
    z:=sqrt(power(x,2) + power(y,2)) - 0.00002 *sin(y * x_pi); 
    theta:= atan2(y, x) - 0.000003 * cos(x * x_pi); 
	  gcj_point:=ST_SetSRID(ST_MakePoint(z *cos(theta),z *sin(theta)),4326); 
	 return gcj_point;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09togcj02_polygon"("source_geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
    target_parts    geometry[];
    source_npoints  integer;
    single_line     geometry;
    single_line_trans geometry;
    single_polygon  geometry;
    final_geom      geometry;

BEGIN
    IF ST_GeometryType(source_geom) != 'ST_Polygon' THEN
        RETURN null;
    END IF;

    FOR single_polygon IN SELECT ST_ExteriorRing ((st_dumprings($1)).geom) as geom LOOP				
        source_npoints := ST_NPoints(single_polygon); 
        single_line  := ST_RemovePoint(single_polygon, source_npoints - 1);  
        single_line_trans := geoc_bd09togcj02_line(single_line);  
        target_parts := array_append(target_parts, ST_AddPoint(single_line_trans, ST_PointN(single_line_trans, 1)));  
    END LOOP;
    SELECT ST_MakePolygon(target_parts[1], target_parts[2:array_upper(target_parts, 1)]) INTO final_geom;  

    RETURN final_geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09towgs84"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
	IF st_srid(geom) != '4490' and  st_srid(geom) != '4326'THEN
        RETURN null;
	end if;
 	CASE ST_GeometryType(geom)
    	when 'ST_LineString' then 
			return geoc_bd09towgs84_line(geom);
		when 'ST_MultiLineString' then 
			return geoc_bd09towgs84_multiline(geom);
		when 'ST_Point' then 
			return geoc_bd09towgs84_point(geom);
		when 'ST_MultiPoint' then 
			return geoc_bd09towgs84_multipoint(geom);
		when 'ST_Polygon' then 
			return geoc_bd09towgs84_polygon(geom);
		when 'ST_MultiPolygon' then
			return geoc_bd09towgs84_multipolygon(geom);
		ELSE
     		RETURN null;
	END CASE;
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09towgs84_line"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		p_p     geometry;
		p_t     geometry;
		z_t     geometry;
		i       int;
BEGIN
    i:=1;
	while i <= st_npoints(geom) LOOP
		p_p := st_pointn(geom,i);
    	p_t := geoc_bd09towgs84_point(p_p);
		geom:=st_setpoint(geom,i-1,p_t);
		i:=i+1;
	end LOOP;
	return geom;

END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09towgs84_multiline"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom $Q$ using geom LOOP
	  	transform_i :=geoc_bd09towgs84_line(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;
CREATE OR REPLACE FUNCTION "public"."geoc_bd09towgs84_multipoint"("geom" "public"."geometry")
  RETURNS "public"."geometry" AS $BODY$
DECLARE
		i                 geometry;
		transform_i       geometry;
		multiArr          geometry[]; 
	
BEGIN
    multiArr:='{}'::geometry[];
	for i in EXECUTE $Q$ select (st_dump($1)).geom $Q$ using geom LOOP
	  	transform_i :=geoc_bd09towgs84_point(i);
		multiArr := array_append(multiArr, transform_i);
	end LOOP;
	return st_multi(ST_Union(multiArr));
END;
$BODY$
  LANGUAGE plpgsql VOLATILE
  COST 100;

使用方法

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
-- 如果转换后结果为null，查看geom的srid是否为4326或者4490
--WGS84转GCJ02
select geoc_wgs84togcj02(geom) from test_table
--GCJ02转WGS84
select geoc_gcj02towgs84(geom) from test_table

--WGS84转BD09
select geoc_wgs84tobd09(geom) from test_table
--BD09转WGS84
select geoc_bd09towgs84(geom) from test_table

--CGCS2000转GCJ02
select geoc_cgcs2000togcj02(geom) from test_table
--GCJ02转CGCS2000
select geoc_gcj02tocgcs2000(geom) from test_table

--CGCS2000转BD09
select geoc_cgcs2000tobd09(geom) from test_table
--BD09转CGCS2000
select geoc_bd09tocgcs2000(geom) from test_table

--GCJ02转BD09
select geoc_gcj02tobd09(geom) from test_table
--BD09转GCJ02
select geoc_bd09togcj02(geom) from test_table
发布：2022-05-22 15:13:33
修改：2023-01-30 21:14:16
链接：https://meethigher.top/blog/2022/gis-science/
标签：gis 
付款码 打赏 分享
访客数calculating人
shift+ctrl+1可控制目录显示