Geo::Transform 作りました
地図と測量の座標変換を手軽にしたい
Geotools (Java) でやたらとややこしいのが座標変換です。まず第一に、必要なクラスを見つけるのがものすごく難しいです。「奥が深い」と思わせるという意味で、バッドノウハウの温床になっています。ということで、座標変換系機能のグッドラッパーを目指す、Geo::Tranform を作ってみました。
実装詳細
| Geo::import_epsg_crs(epsg_code) | 座標系の EPSG コードを与えると、その座標系の Geotools 実装を返します。 |
| Geo::Transform | 座標変換器クラスです |
| Geo::Transform.new(srs_crs, dst_crs) | 座標変換器のコンストラクタです。変換前の座標系、変換後の座標系という順番で引数を渡します。 |
| Geo::Transform#transform(x, y) | 変換前の座標系での座標 (x, y) を変換後の座標系の座標に変換したものを配列として返します。 |
Geotools で3次元の座標の変換はできなかったかどうか、あとで確認したいと思います。（当面3次元の座標変換を扱うつもりはないので、私の作業の中では優先度は低いです。）
動作チェック
ix = Geo::import_epsg_crs(2451) # EPSG:2451 - 平面直角座標系 IX 系 wgs84 = Geo::import_epsg_crs(4326) # EPSG:4326 - WGS84 t = Geo::Transform.new(ix, wgs84) # IX 系から WGS84 への座標変換器 t_inv = Geo::Transform.new(wgs84, ix) # WGS84 から IX 系への座標変換器 pt = t.transform(0, 0) # IX 系の原点を WGS84 に座標変換 pt_inv = t_inv.transform(pt[0], pt[1]) # その点を IX 系に戻す。元に戻るか? print "IX origin is #{pt.inspect} in WGS84\n" print "#{pt_inv.inspect} must be (0, 0)\n"
に対し、
$ ruby geotools.rb IX origin is [139.833333333333, 35.9999999991025] in WGS84 [0.0, 0.0] must be (0, 0) $ jruby geotools.rb IX origin is [139.8333333333333, 35.99999999910247] in WGS84 [0.0, 0.0] must be (0, 0)
という結果を戻してくれます。
この結果は、以前 Ruby + OGR で同様の計算をしたときの結果:
"POINT (139.833333333333428 36.000000000000036)"
http://d.hatena.ne.jp/hfu/20070704/1185663573
と同じと言ってよいと考えています。表示上の（？）細かな違いの評価と、その原因についてはあとで考えたいと思います。
実装
現時点での geotools.rb の実装は以下のようになっています。
# this code is under development and subject to major change. require 'iconv' module Geo # Geo::Tools module, to include nesessary classes from Geotools module Tools QUALIFIED_NAMES = %w{java.lang.String java.lang.Integer java.lang.Double java.lang.Long java.io.File org.geotools.data.shapefile.ShapefileDataStore org.geotools.feature.AttributeTypeFactory org.geotools.feature.FeatureTypeBuilder org.geotools.feature.type.GeometricAttributeType com.vividsolutions.jts.io.WKTReader org.geotools.referencing.crs.EPSGCRSAuthorityFactory org.geotools.referencing.operation.DefaultCoordinateOperationFactory org.geotools.geometry.DirectPosition2D} begin require 'rjb' QUALIFIED_NAMES.each do |qn| sn = qn.split('.').last module_eval "#{sn} = Rjb::import('#{qn}')" end IMPLEMENTATION = 'rjb' rescue LoadError require 'java' QUALIFIED_NAMES.each do |qn| include_class qn end IMPLEMENTATION = 'java' end end # Geo module variables @@wkt_reader = nil @@epsg_crs_authority_factory = nil # Geo module 'good-wrapper' / 'Grossklasstum' classes class Reader ## TODO: implement attribute whitelist filtering (for better performance) def Reader::foreach(shapefile, sjis_workaround = false) if(Tools::IMPLEMENTATION == 'java' && sjis_workaround) raise "sjis_workaround for JRuby is not implemented." end store = Tools::ShapefileDataStore.new(Tools::File.new(shapefile).toURL) iter = store.getFeatureSource.getFeatures.features feat_type = store.getFeatureSource.getSchema attr_names = [] feat_type.getAttributeCount.times do |i| attr_names << feat_type.getAttributeType(i).getName end while(iter.hasNext) feat = iter.next attrs = {} feat.getNumberOfAttributes.times do |i| attr = feat.getAttribute(i) if Tools::IMPLEMENTATION == 'rjb' if attr.getClass.equals(Tools::Integer) attr = attr.intValue elsif attr.getClass.equals(Tools::Double) attr = attr.doubleValue elsif attr.getClass.equals(Tools::String) if sjis_workaround attr = Iconv.conv('UTF-8', 'Shift_JIS', attr.getBytes('iso-8859-1')) else attr = attr.toString end elsif attr.getClass.equals(Tools::Long) attr = attr.longValue end end attrs[attr_names[i]] = attr end attrs.delete('the_geom') yield feat.getDefaultGeometry, attrs end iter.close end end class Writer def Writer::open(shapefile) w = Writer.new(shapefile) yield w w.close end def initialize(shapefile) @shapefile = shapefile @writer = nil @first = true end def setup(geom, attrs) attrs.delete('the_geom') ftb = Tools::FeatureTypeBuilder.newInstance(@shapefile) attrs.each do |key, value| if value.methods.include?('_classname') attr_class = value.getClass elsif value.class == String attr_class = Tools::String elsif value.class == Fixnum attr_class = Tools::Integer elsif value.class == Float attr_class = Tools::Double else raise "attribute #{key} has unrecognizable class #{value.class}" end ftb.addType(Tools::AttributeTypeFactory.newAttributeType(key, attr_class)) end if geom.class == String geom = import_wkt_geometry(geom) end ftb.setDefaultGeometry(Tools::GeometricAttributeType.new('the_geom', geom.getClass, true, nil, nil, nil)) ft = ftb.getFeatureType store = Tools::ShapefileDataStore.new(Tools::File.new(@shapefile).toURL) store.createSchema(ft) @writer = store.getFeatureWriter(@shapefile, store.getFeatureSource(@shapefile).getTransaction) @first = false end private :setup def write(geom, attrs) setup(geom, attrs) if @first feat = @writer.next if geom.class == String geom = import_wkt_geometry(geom) end feat.setDefaultGeometry(geom) attrs.each do |key, value| feat.setAttribute(key, value) end @writer.write end def close @writer.close unless @writer == nil end end class Transform def initialize(src_crs, dst_crs) cof = Geo::Tools::DefaultCoordinateOperationFactory.new @co = cof.createOperation(src_crs, dst_crs) @mt = @co.getMathTransform end def transform(x, y) # z? r = Geo::Tools::DirectPosition2D.new @mt.transform(Geo::Tools::DirectPosition2D.new(x, y), r) return r.x, r.y end ## TODO: def transform(geom) ## TODO: accessor to @mt or @co end # Geo module convenient methods def Geo::import_wkt_geometry(wkt) @@wkt_reader = Geo::Tools::WKTReader.new if @@wkt_reader == nil @@wkt_reader.read(wkt) end def Geo::import_epsg_crs(epsg_code) @@epsg_crs_authority_factory = Geo::Tools::EPSGCRSAuthorityFactory.new if @@epsg_crs_authority_factory == nil if epsg_code.class == Fixnum return @@epsg_crs_authority_factory.createCoordinateReferenceSystem("EPSG:#{epsg_code}") elsif epsg_code.class == String return @@epsg_crs_authority_factory.createCoordinateReferenceSystem(epsg_code) else raise "Geo::import_epsg_crs: can not handle epsg_code = #{epsg_code}" end end def dms2dec(d, m, s) d + m / 60.0 + s / 3600.0 end def dec2dms(dec) #TODO raise "not implemented." end end # ad hoc tests ## TODO: better separate tests as unit tests. ix = Geo::import_epsg_crs(2451) # EPSG:2451 - 平面直角座標系 IX 系 wgs84 = Geo::import_epsg_crs(4326) # EPSG:4326 - WGS84 t = Geo::Transform.new(ix, wgs84) # IX 系から WGS84 への座標変換器 t_inv = Geo::Transform.new(wgs84, ix) # WGS84 から IX 系への座標変換器 pt = t.transform(0, 0) # IX 系の原点を WGS84 に座標変換 pt_inv = t_inv.transform(pt[0], pt[1]) # その点を IX 系に戻す。元に戻るか? print "IX origin is #{pt.inspect} in WGS84\n" print "#{pt_inv.inspect} must be (0, 0)\n"
