描述

輸入

一系列坐標點

算法——CrimeStat

a)? ? ? ?The Y-axis is rotated clockwise through an angle?, where（順時針旋轉）

第四版(有錯誤)

其中

文獻(?Efficacy of Standard Deviational Ellipses in the Application of Criminal Geographic Profiling? )_選用

b)? ? ??Two standard deviations are calculated, one along the transposed X-axis and one

along the transposed Y-axis:

第四版公式(少2)

第三版公式(錯誤)

文獻(Efficacy of Standard Deviational Ellipses in theApplication of Criminal Geographic Profiling? )_選用

where ? and? are the means of X and Y respectively,? is the angle (in radians), and N is the number of points.?

Note, again, that 2 is subtracted from the number of points in both denominators to produce an unbiased estimate of the standard deviational ellipse since there are two constants from which the distance along each axis is measured

c)? ? ?The X-axis and Y-axis of the ellipse are defined by:

d)? ? ?The area of the ellipse is

算法分解

輸出一個橢圓，需有關鍵參數：

圓心、旋轉角度、XY軸長度

輸出

Question&Answer

1. 計算的X,Y是地理坐標還是投影坐標

SQL函數

輸入id 數組 {[? ? ? ? ? ? ? ? ?]}

通過id 獲得 地理位置geography

geography 獲得經緯度? ?

計算中心點X,Y

計算偏角

計算半徑長

畫圓公式

測試

select *from crimestat.v_incidents as PAwhere PA.sid in (select sid from crimestat.v_incidents limit 100 ) --按需獲得數據

select ST_Transform(地理位置::geometry,4326) from crimestat.v_incidents limit 100;——獲得地理坐標

select ST_X(地理位置::geometry)::numeric as x,ST_Y(地理位置::geometry)::numeric as y from crimestat.v_incidents limit 100;——從地理坐標獲得 經(x)緯(y)度?

select ST_X(geom)::numeric as x,ST_Y(geom)::numeric as y? from? ?(select ST_Transform(ST_SetSRID(地理位置::geometry,4326),32651) as geom from crimestat.v_incidents limit 100) as D;

——獲得投影后的公里網坐標x, y.（蘇州東經119°55′～121°20′，北緯30°47′～32°02′之間.https://spatialreference.org/ref/epsg/wgs-84-utm-zone-51n/——EPSG：32651）

--獲得投影后的公里網坐標,輸入參數：投影坐標系EGPS,地理坐標create function getprj( coord int, geo geography)returns table(x numeric,y numeric) as $$select ST_X(geom)::numeric as x,ST_Y(geom)::numeric as yfrom (select ST_Transform(ST_SetSRID($2::geometry,4326),$1) as geom) as D;$$ language sql;select getprj(32651::integer,地理位置::geography)from crimestat.v_incidents as PAwhere PA.sid in (select sid from crimestat.v_incidents limit 100);

try 1

--獲得投影后的公里網坐標,輸入參數：投影坐標系EGPS,地理坐標create or replace function getprj( coord int,out x numeric,out y numeric)returns SETOF record as $$select ST_X(geom)::numeric ,ST_Y(geom)::numeric from (select ST_Transform(ST_SetSRID(地理位置::geometry,4326),$1) as geomfrom crimestat.v_incidents as PAwhere PA.sid in (select sid from crimestat.v_incidents limit 100)) as D;$$ language sql;select * from getprj(32651::integer)

try 2

最終結果

CREATE OR REPLACE FUNCTION crimestat.f_analysis_sde(IN uids json)

? ? RETURNS TABLE(gx numeric, gy numeric)

? ? LANGUAGE 'plpgsql'

? ? VOLATILE

? ? PARALLEL UNSAFE

? ? COST 100? ? ROWS 1000

AS $BODY$declare

N integer;

centerX numeric;

centerY numeric;

angle numeric;

lengthX numeric;

lengthY numeric;

begin

--獲得投影后的公里網坐標 臨時表

create temp table prjxy? on commit drop as

(select ST_X(geom)::numeric as x ,ST_Y(geom)::numeric as y

from ( select ST_Transform(ST_SetSRID(地理位置::geometry,4326),32651) as geom from crimestat.v_incidents as PA

? where PA.uid in (select json_array_elements::character varying from json_array_elements($1))--待查詢記錄

? ) as D) ;

--獲得數量N

select count(*)into N from crimestat.v_incidents as PA

? where PA.uid in (select json_array_elements::character varying from json_array_elements($1));--待查詢記錄

--獲得圓心的X

select sum(x) / N into centerX from prjxy;

--獲得圓心的Y

select sum(y) / N into centerY from prjxy;

--獲得旋轉角度

select atan((sum(power((x-centerX),2))-sum(power((y-centerY),2))+sqrt(power((sum(power((x-centerX),2))-sum(power((y-centerY),2))),2)+4*(power(sum((x-centerX)*(y-centerY)),2)))) / (2*sum((x-centerX)*(y-centerY))))::numeric

into angle from prjxy ;

--獲得X半軸長度

select sqrt(2*sum(power((((x-centerX)*cos(angle))-((y-centerY)*sin(angle))),2))/(N-2))::numeric into lengthX from prjxy;

--獲得y半軸長度

select sqrt(2*sum(power((((x-centerX)*sin(angle))+((y-centerY)*cos(angle))),2))/(N-2))::numeric into lengthY from prjxy;

--構建輪廓點

create temp table geoxy on commit drop as select ST_Transform(ST_SetSRID(ST_Point((centerX+lengthX * cos(radians(seg))*cos(angle)+lengthY*sin(radians(seg))*sin(angle))::numeric , (centerY+lengthY * sin(radians(seg))*cos(angle)-lengthX*cos(radians(seg))*sin(angle))::numeric ),32651),4326) as gpoint

from (select row_number() over() as seg from crimestat.v_incidents limit 360) as Pt;

--獲得360個輪廓點的度數

return query SELECT ST_X(gpoint)::numeric as gx,ST_Y(gpoint)::numeric as gy from geoxy;

end;

$BODY$;