meethigher-Gis坐标系4326与3857及高德百度坐标系转换

已于 2023-06-08 21:33:31 修改
阅读量9.2k 收藏 21
点赞数 5
分类专栏： java 文章标签： geotools gis crs
于 2022-05-22 21:54:54 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/qq_30460361/article/details/124917317
版权

参考

  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

<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。

/**
 * 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百度坐标转高德坐标

进行修改后如下

/**
 * 仅适用于中国的点位坐标系转换
 * 各地图API坐标系统比较与转换;
 * <p>
 * WGS84坐标系：即地球坐标系，国际上通用的坐标系。设备一般包含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 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 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);
    }

    /**
     * 转换为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);
    }

    /**
     * 火星坐标系 (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);
    }


    /**
     * 点位信息
     *
     * @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扩展

create extension postgis

创建方法

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;

使用方法

-- 如果转换后结果为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
言成言成啊
关注
  • 5
    点赞
  • 踩
  • 21
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
墨卡托投影（102100）与地理坐标系（4326）互转工具javascript版本
07-29
墨卡托投影（102100）与地理坐标系（4326）互转工具javascript版本
大型地图切片程序--3857坐标系
01-02
一个基于Web墨卡托投影（3857坐标系）的地图切片程序。不需要地图规则，只需要指定地图左上角和右下角的3857坐标系（原点在中心，向左向下为正）。同时融合arcgis切片的周边地图。切割图像大小不限制。切图效率高，...
GIS中EPSG4326与EPSG3857的区别及转换关系-20220830更新
jadesl的专栏
04-11 7692
​ EPSG4326是基于WGS84椭球的经纬度坐标系（BL）； EPSG3857是基于球体的、web墨卡托投影（伪墨卡托投影）的投影坐标系，范围为纬度85度以下，由于google地图最先使用而成为实施标准，至今，大多互联网地图都使用EPSG3857。该投影是等角投影，适合地图定向及导航，纬度越高，面积变形越大。 GIS坐标系转换(EPSG:4326与EPSG:3857相互转换 ...............
Java【代码分享 04】坐标系说明+WGS84 GCJ02 BD09坐标系转换工具+Java坐标系转换及验证源代码分享（粘贴可用
最新发布
2401_84254364的博客
04-28 994
WGS84 是为 GPS 全球定位系统建立的坐标系统，是世界上第一个统一的地心坐标系，因此也被称为大地坐标系原始坐标系。一般通过 GPS 记录仪记录下来的经纬度，就是基于 WGS84 坐标系的数据。Google 和高德地图定位的的经纬度（国外）都是基于WGS84坐标系的；但是在国内是不允许直接用 WGS84 坐标系标注的，必须经过加密后才能使用。
GIS坐标系：epsg 3857 vs 4326
qq_32649321的博客
05-30 1万+
1.EPSG:4326（WGS84）-地理坐标 EPSG4326是基于WGS84椭球的经纬度坐标系。 EPSG4326 is in degrees - 3D sphere，单位是度。 EPSG4326表示方法： EPSG4326使用参考球面或椭球面上的坐标系统。 特点：利于存储，可读性高。缺点：会导致页面变形。 2.EPSG:3857（WGS84 Web Mercator）-投影坐标 EPSG3857是基于球体的、web墨卡托投影（伪墨卡托投影）的投影坐标系。 在openlayer
地图坐标系转换(4326与3857转换)
vagabond_的博客
12-02 1万+
gis坐标系转换 4326与3857转换
地理信息epsg:4326 和epsg:3857
qq_37891961的博客
03-06 1万+
GIS坐标系小知识 可以看到这两种投影是属于同一坐标系，就是wgs84。需要明确一点的是投影是想用不同的方式展示地图，由于地球是个不规则的椭圆，在做不同地区的地图的时候为了显示其对应的地理特征，采用了不同的投影方式，总的来说坐标系是坐标系，投影是投影（投影的目的是为了展示数据）。 我们该如何选择地图的投影了？看一下两种常用投影的区别就知道了。 ...
ArcGIS Viewer for Flex 4326坐标系与12100坐标系坐标互转
weixin_30323631的博客
05-07 74
1 定义以4326坐标系定位的地图 1 <esri:SpatialReference id="lo" wkid="4326" /> 2 3 <esri:Map crosshairVisible="true" width="100%" height="100%" id="myMap" logoVisible="false"> 4 <...
coordzbzh-GIS坐标转换.rar
04-22
支持单点坐标不同坐标系（空间、平面经纬度、投影坐标系间的转换） 支持中英子午线的设置 可以定义文件格式后批量转换输出 gis数据处理必须工具
数据坐标系转换器（软件）
02-21
提供多种格式和坐标系转换，支持DWG、SHP、MDB等格式，百度、火星、WGS-84等坐标系任意转换；转换快速质量有保证。
坐标系转换工具【支持批量】
12-31
【GIS工具箱】集成了坐标拾取器、坐标系转换工具和随机坐标生成工具；已被300人的公司广泛使用。良心资源，童叟无欺！ 打包出售，更实惠！ https://download.csdn.net/download/jingqia/75970880
天地图4326 切片 Tile shelme
09-11
Arcgis server 切片时。可以选择此文档作为地铁地图切片规则，省去了自己手动添加比例尺
WEBGIS知识-坐标系
03-26
WEBGIS知识-坐标系
GIS基础知识 - 坐标系、投影、EPSG:4326、EPSG:3857
imherer的博客
08-21 8817
最近接手一个GIS项目，需要用到 PostGIS，GeoServer，OpenLayers 等工具组件，遇到一堆地理信息相关的术语名词，在这里做一个总结。 1. 大地测量学 (Geodesy) 大地测量学是一门量测和描绘地球表面的学科，也包括确定地球重力场和海底地形。 1.1 大地水准面 (geoid) 大地水准面是海洋表面在排除风力、潮汐等其它影响，只考虑重力和自转影响下的形状，这个形状延伸过陆地，生成一个密闭的曲面。虽然我们通常说地球是一个球体或者椭球体，但是由于地球引力分布不均（因为密度不同...
EPSG:4326和4490两种常见的坐标系统讲解
HGX-DJK
03-17 2988
在于这两种坐标系统的基准面和使用范围，其中EPSG 4326是国际上广泛使用的坐标系统，而EPSG 4490更专门用于中国地区的地理信息系统应用。在进行地理空间数据处理时，应选择适合你数据所在区域和所需精度的坐标系统。
4326 坐标测量不准确_在程序中准确地测量时间
cuma2369的博客
07-26 262
4326 坐标测量不准确It is quite common to measure the time in programs using APIs like clock() and gettimeofday(). We may also want to measure the time “accurately” for certain purposes, such as measuring a s...
openlayers6【十】EPSG:3857和EPSG:4326区别详解
热门推荐
@必意玲
05-29 3万+
1. 写在前面 在openlayers 6【二】Map地图详解，初始化一个map地图 这篇文章中，其实已经涉及到了这一点，就是为什么这里我们要把中心点的坐标用 fromLonLat()方法进行包裹。fromLonLat() 方法是继承自 ol.proj 这个类。 import { fromLonLat } from "ol/proj"; 在 openlayers 中创建map时候会有个view 属性，该属性下面会存在一个center属性。这个属性就是设置中心位置的坐标。如果不设置就没有中心，也毫无意义。p
gis 地理坐标系转换
09-27
GIS地理坐标系转换是指将不同的地理坐标系之间进行转换的过程。在实际应用中，常常需要将坐标信息从一种坐标系转换为另一种坐标系，以满足不同地图或系统的需求。Gcoord是一个处理地理坐标系的JS库，可以用来修正百度地图、高德地图等互联网地图坐标系不统一的问题。它支持不同地理坐标系之间的转换，如WGS84、GCJ-02、BD09、EPSG4326等。使用gcoord库，可以通过调用相应的方法来进行坐标转换，包括transform(input, from, to)方法和CRS方法。通过这些方法，可以实现经纬度和xy坐标系之间的转换，同时也支持设置参数进行不同类型坐标系的转换。

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值