モジュール:地図/正射図法
表示
< モジュール:地図
モジュールの解説[作成]
return function(size,dp,vc)
local rd=6370000/dp
local u=math.pi/180
local c={{math.cos(vc[1]*u),math.sin(vc[1]*u)},{math.cos(vc[2]*u),math.sin(vc[2]*u)}}
local function invert(zq)
local vp={(zq[1]-size[1]/2)/rd,-(zq[2]-size[2]/2)/rd}
local vq=math.sqrt(1-vp[1]*vp[1]-vp[2]*vp[2])
local y=vq*c[2][1]-vp[2]*c[2][2]
return {math.atan2(vp[1],y)/u+vc[1],math.atan2(vq*c[2][2]+vp[2]*c[2][1],math.sqrt(vp[1]*vp[1]+y*y))/u}
end
local range={invert({0,0})[1],invert({size[1],0})[1],invert({0,size[2]})[2],invert({size[1]/2,0})[2]}
local function convert(vp)
local p={{math.cos(vp[1]*u),math.sin(vp[1]*u)},{math.cos(vp[2]*u),math.sin(vp[2]*u)}}
return {
p[2][1]*(p[1][2]*c[1][1]-p[1][1]*c[1][2])*rd+size[1]/2,
(p[2][1]*(p[1][1]*c[1][1]+p[1][2]*c[1][2])*c[2][2]-p[2][2]*c[2][1])*rd+size[2]/2
}
end
local function clip(inner,ps)
local result=ps
for k=0,1 do
for j=0,1 do
local v=(j==0 and 0-2) or size[k+1]+2
local ps=result
local list={}
for i=1,#ps do
local p={ps[i],ps[i%#ps+1]}
if (p[1][k+1]>=v)==(j==0) then
table.insert(list,p[1])
end
if not (inner and i==#ps) then
if (p[1][k+1]>=v)~=(p[2][k+1]>=v) then
local r={}
r[k%2+1]=v
r[(k+1)%2+1]=p[1][(k+1)%2+1]+(v-p[1][k%2+1])*(p[2][(k+1)%2+1]-p[1][(k+1)%2+1])/(p[2][k%2+1]-p[1][k%2+1])
table.insert(list,r)
end
end
end
result=list
end
end
return result
end
return {
convert=convert,
convertMulti=function(inner,ps)
local t={}
for _,p in ipairs(ps) do
table.insert(t,convert(p))
end
return clip(inner,t)
end,
contains=function(rect)
return rect[1]<range[2] and range[1]<rect[2] and rect[3]<range[4] and range[3]<rect[4]
end,
getSize=function()
return size
end,
getDp=function()
return dp
end
}
end