Convert WGS84 coordinates to GCJ02 coordinates in PostgreSQL

wgs2gcj.sql

create type coord as (x double precision, y double precision);
create type geocoord as (lat double precision, lon double precision);

create or replace function transform(x double precision, y double precision, out lat double precision,
                                     out lon double precision)
  language plpgsql as
$$
declare
  xy   double precision;
  absX double precision;
  xPi  double precision;
  yPi  double precision;
  d    double precision;
begin
  xy = x * y;
  absX = sqrt(abs(x));
  xPi = x * pi();
  yPi = y * pi();
  d = 20.0 * sin(6.0 * xPi) + 20.0 * sin(2.0 * xPi);

  lat = d + 20.0 * sin(yPi) + 40.0 * sin(yPi / 3.0);
  lon = d + 20.0 * sin(xPi) + 40.0 * sin(xPI / 3.0);
  lat = lat + 160.0 * sin(yPi / 12.0) + 320.0 * sin(yPi / 30.0);
  lon = lon + 150.0 * sin(xPi / 12.0) + 300.0 * sin(xPi / 30.0);
  lat = lat * 2.0 / 3.0;
  lon = lon * 2.0 / 3.0;
  lat = lat + -100.0 + 2.0 * x + 3.0 * y + 0.2 * y * y + 0.1 * xy + 0.2 * absX;
  lon = lon + 300.0 + x + 2.0 * y + 0.1 * x * x + 0.1 * xy + 0.1 * absX;
end;
$$;

create or replace function delta(lat double precision, lng double precision, out dLat double precision,
                                 out dLon double precision)
  language plpgsql as
$$
declare
  earthR    double precision;
  ee        double precision;
  radLat    double precision;
  magic     double precision;
  sqrtMagic double precision;
  tcoord    geocoord;
begin
  earthR = 6378137.0;
  ee = 0.00669342162296594323;
  tcoord = transform(lng - 105.0, lat - 35.0);
  dLat = tcoord.lat;
  dLon = tcoord.lon;
  radLat = lat / 180.0 * pi();
  magic = sin(radLat);
  magic = 1 - ee * magic * magic;
  sqrtMagic = sqrt(magic);
  dLat = (dLat * 180.0) / ((earthR * (1 - ee)) / (magic * sqrtMagic) * pi());
  dLon = (dLon * 180.0) / (earthR / sqrtMagic * cos(radLat) * pi());
end;
$$;

create or replace function wgs2gcj(lat_wgs double precision, lon_wgs double precision, out lat_gcj double precision,
                                   out lon_gcj double precision)
  LANGUAGE plpgsql as
$$
declare
  d geocoord;
begin
  d = delta(lat_wgs, lon_wgs);
  lat_gcj = lat_wgs + d.lat;
  lon_gcj = lon_wgs + d.lon;
end;
$$;

select wgs2gcj(39.92229, 116.34373);