【SU Ruby教程】几何与变换(4):地理数据

Apiglio

共 13772字,需浏览 28分钟

 · 2021-03-07


在了解了普遍性的空间数据表达以及微观的空间关系判断与变换之后,本篇作为有关 Geom 模块的最后一篇,需要介绍一种专门针对宏观地理空间数据的表达方式。涉及了 Geom:: LatLong 和 Geom:: UTM 两个类,分别用来表示地理坐标系的经纬度信息和投影坐标系的水平坐标。


这两个类的功能是将这类地理信息换算成模型中的坐标,而这个换算过程需要当前模型已经设置了地理信息之后才能有效使用,因此本篇在介绍两个地理数据类之间,需要先介绍如何设置地理位置以及 SketchUp 模型如何表示地理位置。




几何与变换(4):地理数据



【本期目录】

(1)当前模型的地理位置

①原点地理坐标

②设置地理位置的方法

③地理位置的存储形式

(3)UTM类

何谓UTM

Geom::UTM

③涉及UTM的转换

(2)经纬度类

①Geom::LatLong

②坐标转为经纬度

③经纬度转为坐标

(4)投影坐标系补遗

①投影坐标的限度

②赤道地区的困境

③自定义坐标转换


(1)当前模型的地理位置


①原点地理坐标


SketchUp模型中关联模型坐标和地理空间,是通过记录原点的地理坐标实现的。在早期版本中还有坐标北方的设置,新版本统一以y轴(绿轴)正方向为北方。


②设置地理位置的方法


有两种方法设置SU模型的地理位置,一种是通过软件界面进行设置,另一种则是直接通过ruby代码来修改地理位置。


(i)通过软件界面设置


第一步:在菜单中选择模型信息,打开模型信息的设置窗口:

第二步:在窗口的地理位置选项卡中选择添加位置:

点击后就可以在地图中选择模型的范围:

第三步:选择具体范围,定位时不接受过大的影像图片参考,不过可以在定位之后继续添加更多影像。

选择确认之后,这张影像图就会作为图片参考以锁定的状态插入到模型中,作为位置参照。同时,模型信息中的位置信息也一并更新。

另外,也可以在第二步中选择“手动设置位置”:

手动设置位置则不会有参考影像图辅助定位,但是定位信息已经储存在当前模型中了。


(ii)通过ruby代码设置


使用代码可以直接定义坐标原点的地理信息:

Sketchup.active_model.shadow_info["Latitude"]=26.0Sketchup.active_model.shadow_info["Longitude"]=119.0Sketchup.active_model.shadow_info["Country"]="China"Sketchup.active_model.shadow_info["City"]="Foochow"

代码的设置方法与“手动设置位置”的窗口设置并无不同,不过需注意经纬度中,以东经为正、西经为负、北纬为正、南纬为负。国家与城市的信息不影响坐标定位。


③地理位置的储存形式


根据上文第二种设置地理位置的方法,就可以了解到 Sketchup. active_model. shadow_info 就可以访问到当前模型的地理位置,因此这个方法返回的 Sketchup:: ShadowInfo 类实例就是可以理解为当前模型地理信息的存储位置。


如果将阴影设置转为哈希类型就能看到以下内容:

Sketchup.active_model.shadow_info.to_h#以下为结果{ "City"=>"Foochow", "Country"=>"China", "Dark"=>45, "DayOfYear"=>183, "DaylightSavings"=>false, "DisplayNorth"=>false, "DisplayOnAllFaces"=>true, "DisplayOnGroundPlane"=>true, "DisplayShadows"=>true, "EdgesCastShadows"=>false, "Latitude"=>26.0, "Light"=>80, "Longitude"=>119.0, "NorthAngle"=>0.0, "ShadowTime"=>'2021-07-02 19:07:54 +0800', "ShadowTime_time_t"=>1625224074, "SunDirection"=>'Vector3d(0.239016, -0.0380432, 0.97027)', "SunRise"=>'2021-07-02 13:19:56 +0800', "SunRise_time_t"=>1625203196, "SunSet"=>'2021-07-03 02:55:53 +0800', "SunSet_time_t"=>1625252153, "TZOffset"=>8.0, "UseSunForAllShading"=>false}#其中单引号包含的内容为输出结果,不是可识别的hash格式。


由于 Sketchup:: ShadowInfo 类最初设计的用途是作为计算阴影的依据,因此大部分字段均是与阴影计算有关的设置,只不过在计算太阳高度角和日照时间时需要用到经纬度信息,就一并记录下来了。除了原点经纬度坐标以外其他的设置会在之后有关文档设置的篇幅中逐一罗列。


除了 ShadowInfo 类以外,还有一个存储模型坐标位置,称为 “GeoReference”。这个设置用于记录地理影像图的定位。


这个 GeoReference 选项属于 Sketchup:: AttributeDictionary 类,其容器类为 Sketchup:: AttributeDictionaries 类,是之后图元属性篇的重要内容,而在此处只需要大概了解它能够存储什么数据即可。


使用与 ShadowInfo 类相似的查看方式可以看到,在不同几种设置地理位置的方法之后,这个设置有不同的表现。

Sketchup.active_model.attribute_dictionaries["GeoReference"].to_h#没有设置定位时的结果{"GeoReferenceNorthAngle"=>0.1557307335588798, "Latitude"=>40.018309, "LocationSource"=>"Manual", "Longitude"=>-105.242139, "ModelTranslationX"=>-18871519.616960514, "ModelTranslationY"=>-174402260.58333763, "UsesGeoReferencing"=>false}#手动设置定位后的结果{"GeoReferenceNorthAngle"=>359.35706222732586, "Latitude"=>40.0, "LocationSource"=>"Manual", "Longitude"=>118.0, "ModelTranslationX"=>-23045687.47413295, "ModelTranslationY"=>-174340002.5487961, "UsesGeoReferencing"=>true}#导入影像图定位后的结果{"GeoReferenceNorthAngle"=>358.7176901069756, "Latitude"=>29.597285947065465, "LocationSource"=>"Google Earth", "Longitude"=>119.59426497463414, "ModelHereState"=>"", "ModelHereZoom"=>15, "ModelTranslationX"=>-29577134.01149756, "ModelTranslationY"=>-129014701.47825345, "ModelTranslationZ"=>-2542.457294041363, "TimeStamp"=>1614759574, "UsesGeoReferencing"=>true, "ZValueCentered"=>-2542.457294041363}


当完全没有设置地理位置信息时: "LocationSource" 字段为 "Manual",表示没有软件给定的影像图; "UsesGeoReferencing" 字段为 false,表示没有启用此选项的经纬度。未启用时, GeoReference 设置的经纬度指向 SketchUp 在Boulder的公司地址。此时模型的地址由 ShadowInfo 确定,通常中文版的默认位置是北京故宫博物院附近。


如果手动设置了地理信息, "UsesGeoReferencing" 字段则为 true,经纬度信息与 ShadowInfo 一致;如果导入了SU提供的影像底图, "LocationSource" 字段为会对底图来源进行说明。


使用软件界面修改模型地理位置时,会同时修改 ShadowInfo 和 GeoReference 两个设置。如果通过代码分别修改,造成两个设置不同时,以 ShadowInfo 的经纬度为转换依据。


(2)经纬度类


①Geom:: LatLong


类似于 Point3d,  Geom 模块定义有 LatLong 类专门储存经纬度信息,尽管这个类型完全可以用 [long, lat] 来替换。

pos = Geom::LatLong.new(34.5108.5)pos.latitude#>> 34.5pos.longitude#>> 108.5puts pos#>> LatLong(34.50000北, 108.50000东)puts pos.to_utm#>> UTM(49 S 270465.57254 3820434.69231)pos.to_a#>> [34.5, 108.5]


构造函数同样接受以下两种参数:

pos1 = Geom::LatLong.new([34.5108.5])pos2 = Geom::LatLong.new(pos)


与 Point3d 不同, LatLong 没有诸如 .distance 这样的方法需要调用,唯一不能用数组类型代替的是 .to_utm 方法,但是 UTM 类本身使用范围并不广,因此 LatLong 类的使用不像 Point3d 那样有必要,以至于转换坐标的方法都不将此类的实例作为可接受的参数类型。


②坐标转为经纬度


因为模型有地理位置信息,所以坐标和经纬度之间转换的方法都是在Sketchup:: Model类中。

pt=Geom::Point3d.new([0,0,0])Sketchup.active_model.point_to_latlong(pt).to_a#>> [119.00000000013267, 25.999999999986425, 0.0]


Model 类中的 .point_to_latlong 方法需要一个 Geom:: Point3d 或 Array 类,而它的输出结果并不是Geom:: LatLong,而是 Point3d。使用点坐标中的x作为经度,y则作为纬度。此方法在SketchUp 6.0之前返回的确实是 Geom:: LatLong,不知出于什么目的修改成了这样。


如果需要返回 LatLong 类,可以自己定义一个:

module Sketchup  class Model    def point2latlong(pt)      return Geom::LatLong.new(self.point_to_latlong(pt).to_a[0..1].reverse)    end  endend


③经纬度转为坐标


同理,经纬度转为坐标也是通过 Sketchup:: Model 的方法:

Sketchup.active_model.latlong_to_point([119.2, 29.3])#>> (13505.20874m, 366016.259858m, 0m)


它同样不接收 LatLong 类的实例,而是需要 [long, lat] 数组。这一点需要专门说明, LatLong 类的 .to_a 方法返回的是 [lat, long] 格式的数组,这意味着如果要使用 LatLong 类,则需要使用 .to_a .reverse 的表述。


同样可以自己定义一个接受 LatLong 实例的版本:

module Sketchup  class Model    def latlong2point(latlong)      self.latlong_to_point(latlong.to_a.reverse)    end  endend


(3)UTM类


①何谓UTM


UTM,即通用横轴墨卡托格网系统(Universal Transverse Mercartor Grid System),与国内地形图常用的高斯-克吕格投影都是等角横轴圆柱投影的格网系统,区别在于前者是割投影,后者是切投影。GPS 的坐标系统为 WGS-84,而 WGS-84 最常见的投影方式就是UTM。


UTM将地球沿经线划分成60个分带,每个分带经度跨度为6°。然后将分带按照一定的投影规则投影到平面上,这个投影规则保证投影后的空间不会发生角度上的失真,因此是等角投影。

(这里展示了投影方式,但是忽略了投影面与地球表面相割的细节。所以准确的说,这个示意图是高斯克-吕格投影而非UTM。


根据 wiki.gis.com 中的描述,UTM的60个分度带从纬度180°向西开始编号,180°-174°W 为1,174°W-168°W 为2,以此类推。根据纬度可以进一步将分度带划分成22个区域,两极所在的分区纬度跨度为10°,其余分区高度为8°。从南纬84°到北纬84°分别用字母CX命名,中间跳过IO两个易与数字混淆的字母。这些字母中MN之间为赤道,可以很方便的记为“N字母开始为Northern Hemisphere”。因此通过经纬网可以将地表划分成数字加字母的格网,例如 120°E-126°E,24°N-32°N 的范围编号为51R。特殊地,纬线圈80°以内的圆形区域根据东西半球划分成两个区,南极为AB区,北极为YZ区。X带的上界从 80°N 上移到 84°N 以使YZ两个区只包含海域。高纬度地区出于更好地表现海陆轮廓的考虑,还有几处例外,这里不展开。


(图片来源:http://wiki.gis.com/wiki/index.php/File:Utm-zones.jpg)


不过 Sketchup 中的 UTM 似乎并不是按照这样的规则来划分的。总体的分区思想不变,但是取消了X区上限北移的做法,而是删除了Z区,在 80°N-84°N 之间增加60个Y分区。同时根据海陆轮廓进行的其他例外的调整和并区的作法也没有使用。


除了地图分幅,UTM 更重要的作用是经纬度转平面坐标,坐标以米为单位,这样才能够使经纬度坐标之间能够计算距离。一个6°分度带上,以中央经线和赤道的交点为实际坐标原点,但是为了使投影坐标都是正数,进行原点偏移。由于赤道上3°的长度为333.96km左右,所以假定原点向西偏移500km可以保证所有横坐标都为正数。如果是南半球的投影,则将原点继续先南偏移10000km,这样也可以使除了A区和B区以外的南半球所有纵坐标都为正数。极地分幅有特殊的规则,其中极点坐标为 (2000000, 2000000),也是出于同样的理由,在这里不展开。


②Geom:: UTM


所以 Geom 模块中的 UTM 类,类似于 LatLong 类,记录地理空间中一点的坐标。而这个点坐标有四个数据,分别是带号、格网图号的字母、平面直角坐标的两个值。


例如以下定义了一个北半球51带的 (500000, 0) 点,也就是该带的实际原点,具体来说就是 123°E 0°N。

u = Geom::UTM.new(51,"N",500000.0,0.0)puts u.zone_number  #>> 51puts u.zone_letter  #>> Nputs u.x            #>> 500000.0puts u.y            #>> 0.0puts u.to_latlong#>> LatLong(0.00000北, 123.00000东)


构造方法同样接受以下两种形式的参数:

u1 = Geom::UTM.new([51,"N",500000.0,0.0])u2 = Geom::UTM.new(u)


其中,当 y>1000000 时,utm表示的点在南半球,这与官方文档中的说法似乎有一些冲突:

Note: Valid ranges for #zone_number and #zone_letter are 1-60 and C-X (omitting I and O). Valid ranges for #x and #y are 100000-899999.

http://ruby.sketchup.com/Geom/UTM.html

因此,如果在使用UTM类的时候发现一些不寻常的特性,应当从理解UTM坐标本身入手进行方法上的修改。


③涉及UTM的转换


UTM 涉及与 LatLong 的转换,也涉及与模型坐标的转换。


UTM 与 LatLong 的转换如下:

u = Geom::UTM.new(51,"N",500000.0,0.0)puts u.to_latlong#>> LatLong(0.00000北, 123.00000东)ll = Geom::LatLong.new(34.5108.5)puts ll.to_utm#>> UTM(49 S 270465.57254 3820434.69231)


虽然 .to_latlong 方法主要以直角坐标和带号为依据,但并不能完全忽视 zone_letter 值错误时的影响。

utm1 = Geom::UTM.new(51,"N",500000.0,0.0)utm2 = Geom::UTM.new(51,"X",500000.0,0.0)utm3 = Geom::UTM.new(51,"F",500000.0,0.0)puts utm1.to_latlongputs utm2.to_latlongputs utm3.to_latlong#>> LatLong(0.00000北, 123.00000东)#>> LatLong(0.00000北, 123.00000东)#>> LatLong(90.01823北, 123.00000东)

因此至少要保证 zone_letter 与坐标实际表示的点在赤道同侧,否则就会出现如上 utm3 的超值错误。而这个错误来自于 to_latlong 方法需要通过 zone_letter 来确定是否需要将纵坐标还原回10000km以南,以反映南半球坐标的实际位置。南半球的坐标输入北半球的字母代号也会有相同的问题。


UTM 与模型点坐标的转换同样需要 Sketchup:: Model类:

model = Sketchup.active_modelutm_1 = model.point_to_utm([0,0,0])#>> UTM(50 R 694820.98408 3209634.52687)puts utm_1.to_latlong#>> LatLong(29.00000北, 119.00000东)# 这是模型地理位置信息设置中的经纬度utm_2 = model.point_to_utm([100.m,100.m,0])#>> UTM(50 R 694906.58036 3209749.83908)pt1 = utm_1.to_a[2..3]pt2 = utm_2.to_a[2..3]r=Math.sqrt((pt1[0]-pt2[0])**2+(pt1[1]-pt2[1])**2)i=Math.sqrt(100**2*2)puts r/i#>> r = 143.60929016770177#>> i = 141.4213562373095#>> 1.0154710291896851


从例子中可以看到模型中的坐标点 [0,0,0] 和 [100.m,100.m,0] 转为UTM坐标之后,坐标距离和实际距离相比有一些扩大。这是因为任何一种地图投影都无法同时在平面上兼顾角度变形与面积变形。因而, UTM作为等角投影,在不同纬度上面积变形的程度是不同的。当上述例子中使用与北方成45°角的线段长度来测试变形时,这个变形就十分明显。如果此时将模型的地理位置的纬度进行修改,最后计算出的比值就会发生改变。例如改成 119°E 1°N 时,最终的结果为0.99687。这也能体现 UTM 是割圆柱投影,在靠近赤道的地区由于投影面在地表以下,因此投影面积要比实际面积小。如果是高斯克吕格投影,由于是切圆柱投影,面积只会扩大而不会缩小。


因此SketchUp在处理过大的地理数据时会出现严重的形变问题,因此一定要注意模型的尺寸问题。


可以使用 .point_to_utm 方法将模型点坐标转换为UTM,也同样可以通过 .utm_to_point 方法进行逆向操作:

model = Sketchup.active_modelutm_p = Geom::UTM.new(51,"N",500000.0,0.0)point = model.utm_to_point(utm_p)#>> (-3228265233.004935mm, 445277963.173094mm, 0mm)


(4)投影坐标系补遗


关于投影坐标系还要再简要补充一些内容,关于地理坐标与模型精度的问题。


①投影坐标的限度


SketchUp模型以局部建模见长,并不太适用于大面积的地理要素建模。尤其是SU的平面坐标系以UTM坐标系为基础,如果在较大范围内使用经纬度坐标就会出现严重的形变。


尝试以下代码:

model=Sketchup.active_modelp=[]max_lng=117+3min_lng=117-3max_lat=80min_lat=0for lat in min_lat..max_lat do  p<<(min_lng..max_lng).to_a.collect{|lng|    Geom::LatLong.new(lat,lng)  }endpt=p.collect{|arr|  arr.collect{|latlong|    model.latlong_to_point(latlong.to_a.reverse)  }}for i in 0..(max_lat-min_lat) do  for j in 0..(max_lng-min_lng-1) do    model.entities.add_line(pt[i][j],pt[i][j+1])  endendfor i in 0..(max_lat-min_lat-1) do  for j in 0..(max_lng-min_lng) do    model.entities.add_line(pt[i][j],pt[i+1][j])  endend


此代码用于生成当前地理位置投影的经纬网,范围是 114°E-120°E 0°-80°N。可以得到以下结果:


可以看到在6°范围内,形变没有很明显。但是如果继续扩大经度范围,使得坐标超出UTM分度带。此时仍然可以投影坐标,但是形变程度会逐渐扩大。例如扩大到18°可以得到以下结果:


因此,如果需要使用SketchUp制作更大尺度的三维作品,就需要自定义的投影坐标,这个会在最后提及。


②赤道地区的困境


在讲自定义投影规则之前,需要再补充一个和UTM坐标有关的问题。如果使用上述代码创建完整的一个分度带经纬网,就会导致一个新的错误:

这是因为 .latlong_to_point 方法和 Geom:: UTM 的 .to_latlong 方法一样,存在南半球坐标的转换问题。如果模型的地理范围恰巧地跨赤道两侧,而又需要使用经纬度数据,这无疑是一个灾难,这时候就需要特别地将南半球坐标平移10000km,代码如下:

def latlong2point(latlong)  if latlong.is_a?(Array) then    lat=latlong[1]    lng=latlong[0]  elsif latlong.is_a?(Geom::LatLong) then    lat=latlong.latitude    lng=latlong.longitude  else    return(nil)  end  pt=Sketchup.active_model.latlong_to_point([lng,lat])  if lat<0 then pt[1]-=10000.km end  return ptend
model=Sketchup.active_modelp=[]max_lng=117+3min_lng=117-3max_lat=89min_lat=-89for lat in min_lat..max_lat do p<<(min_lng..max_lng).to_a.collect{|lng| Geom::LatLong.new(lat,lng) }endpt=p.collect{|arr| arr.collect{|latlong|    latlong2point(latlong) }}for i in 0..(max_lat-min_lat) do for j in 0..(max_lng-min_lng-1) do model.entities.add_line(pt[i][j],pt[i][j+1]) endendfor i in 0..(max_lat-min_lat-1) do for j in 0..(max_lng-min_lng) do model.entities.add_line(pt[i][j],pt[i+1][j]) endend

效果如下:



③自定义坐标转换


如果模型需要更大的空间范围,就不能用UTM这种分带投影了,这时如果还需要使用SketchUp,那么就需要考虑自己设置投影规则,我们可以参考 Geom:: UTM 类,创建一个自定义的坐标类,从而实现地理坐标的自定义投影。以下展示一个简单的例子:


暂且称它为“Apiglio墨卡托投影”,其实就是垂直等距的绘制经线和纬线,每个角度在投影面上的长度为1km。


首先在 Geom 模块中追加一个这个类的定义,就像 UTM 类一样,有对应的构造函数和访问内部坐标的方法,也包括了转为 LatLong 类的方法。其中的 .initialize 方法是专门用于声明的,访问这个方法需要使用类名+“.new”。

module Geom  class ApiglioMercator    attr_accessor :lat, :lng    def initialize(_lat,_lng)      @lat,@lng=_lat,_lng    end    def to_a      [@lng,@lat]    end    def to_latlong      Geom::LatLong.new(@lat,@lng)    end  endend


之后继续在 Model 类里头追加转换方法。尽管这个方法并不需要读取 Model 类的 ShadowInfo 设置,但是出于统一,转为模型点坐标的方法都应该放在这个类中:

module Sketchup  class Model    def apigliomercator_to_point(am)      unless am.is_a?(Geom::ApiglioMercator) then return nil end      return Geom::Point3d.new(am.lng*1000.m,am.lat*1000.m,0)    end    def point_to_apigliomercator(point)      pt=Geom::Point3d.new(point)      return Geom::ApiglioMercator.new(pt[1]/1000.0.m,pt[0]/1000.0.m)    end  endend


最后为了展示效果,这里选择了这样一个文件数据,每一行四个浮点数分别表示边线的起讫点经纬度坐标:


然后执行以下代码:

filename="K:/..."#给一个具体的文件地址f=File.readlines(filename)f.collect!{|i|i.split(",")}f.collect!{|edg|edg.collect(&:to_f).each_slice(2).to_a}f.collect!{|edg|  edg.collect{|pt|    Sketchup.active_model.apigliomercator_to_point(      Geom::ApiglioMercator.new(pt[1],pt[0])    )  }}f.each{|edg|  Sketchup.active_model.entities.add_line(edg[0],edg[1])}


执行可能需要一些时间,执行完成后得到以下效果:


以上就是本篇教程的全部内容。




最后,如果没有启用 GeoReference,其经纬度指向SketchUp的公司地址,顺着经纬度可以找到这样的位置:


然后对照现在SketchUp的东家——Trimble公司的官网,就可以发现其中一个分部就位于这里:





本文编号:SU-R09

浏览 99
点赞
评论
收藏
分享

手机扫一扫分享

举报
评论
图片
表情
推荐
点赞
评论
收藏
分享

手机扫一扫分享

举报