plsql根据3度投影计算面积

postgis的经纬度数据，为了稍微精确计算面积，通过获取中心点的经度坐标，计算3都投影的epsg，并计算面积（cgcs2000坐标系）。

原来用select ST_AREA(ST_Transform(ST_SetSRID(geom,4326), epsg)) into area;

会非常慢；

更改为直接处理，如下所示：

create or replace function area_3dgree_new(geom geometry)
returns double precision as
$$
declare
area double precision=0;
lon double precision=ST_X(ST_Centroid(geom));
dh integer=floor((lon+1.5)/3);
epsg integer=(dh - 25)+4534;
srid integer=ST_Srid(geom);
begin
	if position('polygon' in lower(GeometryType(geom))) > 0 then
		if (srid<=0 and lon < 360) or srid=4326 then
			srid=4490;
			geom=ST_SetSRID(geom,srid);
		end if;

		if srid in (4490) then
			area= ST_AREA(ST_Transform(geom, epsg));
		else
			area=ST_AREA(geom);
		end if;
	end if;
	return area;
end
$$
language 'plpgsql' stable;

获取2000坐标系三度带带号：

create or replace function epsg_3dgree(geom geometry,onlyGeography bool = true)
returns integer as
$$
declare
srid integer=ST_Srid(geom);
area double precision=0;
lon double precision;
dh integer;
epsg integer;
begin	
	if onlyGeography then
		if not srid in (4326,4490) then
			return srid;
		end if;
	end if;
	if srid in (4326,4490) then
		lon = ST_X(ST_Centroid(geom));
	else
		lon=ST_X(ST_Centroid(ST_Transform(geom,4490)));
	end if;
	dh=floor((lon+1.5)/3);
	epsg=(dh - 25)+4534;
	return epsg;
end
$$
language 'plpgsql' stable;

转换投影：

create or replace function epsg_transform(geom geometry,target_srid integer)
returns geometry as
$$
declare
srid integer=ST_Srid(geom);
begin	
	if target_srid=srid then
		return geom;
	end if;
	if srid=0 or target_srid=0 then
		return ST_SetSRID(geom,target_srid);
	end if;
	if srid=4326 then
		if target_srid=4490 then
			return ST_SetSRID(geom,target_srid);
		else
			return ST_Transform(ST_SetSRID(geom,4490),target_srid);
		end if;
	end if;
	if target_srid=4326 then
		if srid=4490 then
			return ST_SetSRID(geom,target_srid);
		else
			return ST_SetSRID(ST_Transform(geom,4490),target_srid);
		end if;
	end if;
	return ST_Transform(geom,target_srid);
end
$$
language 'plpgsql' stable;