Test used Spatialite transform SRID (BaiDemon function) 使用Spatialite 轉換 SRID 與 我的測試方法

目的:

找出Table 的 SRID 轉換方式 使用Spatialite (舉例:4326 wgs84 轉 tw97/ TM2 zone 121)

測試方式:

使用QGIS 自動轉換座標的功能,例如說我今天匯了兩個不同座標系統的幾何資料,而我QGIS的主要座標系統(桌布的座標系統) 設定為WGS84 的話,只要我會進來的幾何資料不管是TW97或TW67或是其他更多種不管,都會自動轉換為我的桌布主要座標系統,所以會自動的套疊起來,

而今天我們要做兩個檔案,來作測試

Make This Files:

1.Point3826

2.Polygon4326

以上因為QGIS會自動轉換座標,在實作時我套疊再一起,我畫一張圖就看得出來了

　　我們為什麼要這樣做呢,還有我的框是畫在點的邊緣,為什麼要這樣做,我只是要測試轉換後是不是會有誤差,而怎麼測試呢? spatialite就不會像QGIS這樣的人性化,所以不同的座標系系統在作幾何運算時,不會自動轉換,所以我們要將不同的座標系統,轉換成同一個座標系統,再來作運算,轉換完後再來比對說,到底有沒有失真,意思就是出現誤差

所有的步驟如下::

1.作出上圖的兩個圖層,不同的SRID,我是拿 4326 和3826 當作範例

2.使用Spatialite你可以把4326轉成3826 ,都可以拉

3.轉換完後使用幾何運算,within()這個function來看說到底,有沒有和QGIS的套碟內容一樣

工具:

1.QGIS

2.Spatialite

3.Spatialite_Gui

從這裡開始操作步驟::

請先用spatialite 作一個table出來

1.可以使用win command 模式 cd 到spatialite檔案的資料夾下

2.下指令 spatialite testdb.sqlite

3.這時候會自動進入testdb.sqlite , 請下 create table 123 ; (這個是錯誤的指令,但是會產生這個資料庫的實體,我們暫時只需要空的實體資料庫)

4.圖2.請使用Spatialite_gui 連結資料庫後,load 用QGIS作出來的兩個不同SRID的幾何資料

新增說明文字

5. ALTER TABLE testid3826 add column wgs84; (先在這個欄位上加入wgs84的這個column)

6.update testid3826 set wgs84 = transform(geometry,4326) ; (把轉換過的座標寫入wgs84column內)

7.輸入圖內的指令,你可以看到轉換前和轉換後,而介紹一下astext()這個函數就是顯示座標,你也可以打select * from table 看所有,就會知道為什麼要用astext這個function

8.實作幾何運算,使用within()function ,

   8.1請先試試看沒有轉換過的

    select a.* from testid3826 as a

    join testid4326 as b 

    where within(a.geometry,b.geometry)

    答案是,null 就是比對不到

   8.2試試看轉換過的

   select a.* from testid3826 as a

   join testid4326 as b 

   where within(a.wgs84,b.geometry)

   答案正確,你可以自己比對看看喔

以下參考

(以下有一行比對所有srid 的table 是 spatial_ref_sys)

可能有版本上的差異

再比較新版的 你可以看一下 select * from spatial_ref_sys

順便記錄一下,如果 load須要給SRID的話 ,但是不給 給0的意思是::

From 

http://www.gaia-gis.it/spatialite-1.0a/SpatiaLite-tutorial.html

Translate Taiwanese

5. Managing Coordinate Reference Systems and Coordinate Transformation

Each coordinate value, to be fully qualified, needs to be explicitly assigned to some SRID, or Spatial Reference Identifier.

This value identifies the geometry's associated Spatial Reference System that describes the coordinate space in which the geometry object is defined.
As a general rule, different geometries can interoperate in a meaningful way only if their coordinates are expressed in the same Coordinate Reference System.
Very often some kind of coordinate reprojection is required in order to convert different GIS data in an unique, homogeneous, Coordinate Reference System and thus allowing interoperation and integration.

When you try to use two GIS datasets belonging to different Coordinate Reference Systems, entities falls very far because numeric values for coordinates are obviously in two different spaces.

But if you apply some opportune coordinate reprojection for one dataset, so to put it in the same Coordinate Reference Systems of the other one, entities will correctly overlap as expected.

The European Petroleum Survey Group [EPSG] maintains and distributes a large dataset of geodetic parameters describing quite any Coordinate Reference System and Coordinate Transformation used worldwide for GIS data.
SpatiaLite on its own implements SRIDs for any kind of geometric class, and supports the EPSG dataset to identify Coordinate Reference Systems.
The test.db already contains the epsg table, so we can quickly start an SQLite session:

sqlite> SELECT * FROM epsg LIMIT 5;

srid	name	proj_params
2000	Anguilla 1957 / British West Indies Grid	+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs
2001	Antigua 1943 / British West Indies Grid	+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +units=m +no_defs
2002	Dominica 1945 / British West Indies Grid	+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=725,685,536,0,0,0,0 +units=m +no_defs
2003	Grenada 1953 / British West Indies Grid	+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=72,213.7,93,0,0,0,0 +units=m +no_defs
2004	Montserrat 1958 / British West Indies Grid	+proj=tmerc +lat_0=0 +lon_0=-62 +k=0.9995000000000001 +x_0=400000 +y_0=0 +ellps=clrk80 +towgs84=174,359,365,0,0,0,0 +units=m +no_defs

sqlite>

1. each Coordinate Reference System is uniquely identified by a srid value.
2. and is qualified by an intelligible name
3. and has a lot of [quite obscure] geodetic projection parameters
   Don't worry about this ... you can just get them as black magic for now ...
   If you want to learn more about them, you can usefully read the followings: http://en.wikipedia.org/wiki/Geodetic and ftp://ftp.remotesensing.org/proj/OF90-284.pdf

sqlite> SELECT DISTINCT Srid(GaiaGeometry) FROM Towns;

Srid(GaiaGeometry)
32632

sqlite>

1. the SpatiaLite SRID() function allows you to identify the srid value that identifies any kind of geometry.
2. all geometries in the test.db have a srid value of 32632.
   You can check this by yourself.

sqlite> SELECT DISTINCT SRID(Towns.GaiaGeometry), epsg.name
     ...> FROM Towns, epsg
     ...> WHERE SRID(Towns.GaiaGeometry) = epsg.srid;

SRID(Towns.GaiaGeometry)	epsg.name
32632	WGS 84 / UTM zone 32N

sqlite>

1. you can perform a simple join to discover what srid 32632 really means.
2. all geometries in the test.db pertains to the UTM zone 32N Coordinate Reference System.

Now suppose that for any good reason you absolutely need to transform coordinates, in such a way to obtain a new dataset pertaining to the WGS 84 Coordinate Reference System.
You can do this quite complex task by simply using a plain SQL query.

sqlite> BEGIN;
sqlite> ALTER TABLE Towns ADD COLUMN wgs84 BLOB;
sqlite> UPDATE Towns SET wgs84 = Transform(GaiaGeometry, 4326);
sqlite> COMMIT;
sqlite> SELECT AsText(GaiaGeometry), Srid(GaiaGeometry),
     ... > AsText(wgs84), Srid(wgs84) FROM Towns LIMIT 5;

AsText(GaiaGeometry)	Srid(GaiaGeometry)	AsText(wgs84)	Srid(wgs84)
POINT(427002.77 4996361.33)	32632	POINT(8.071929 45.116952)	4326
POINT(367470.48 4962414.5)	32632	POINT(7.324219 44.802838)	4326
POINT(390084.12 5025551.73)	32632	POINT(7.596214 45.374854)	4326
POINT(425246.99 5000248.3)	32632	POINT(8.049028 45.151753)	4326
POINT(426418.89 4957737.37)	32632	POINT(8.070136 44.769244)	4326

sqlite>

1. the srid that identifies WGS 84 is 4326.
   You can check this by yourself.
2. the SpatiaLite Transform() function obtains a new geometry from an original one, applying a coordinate reprojection as required.
3. Now you have two alternative geometries for the Towns table:
   • the original GaiaGeometry column contains geometries in the 32623 - UTM zone 32N Coordinate Reference System.
   • the new wgs84 column contains geometries in the 4326 - WGS84 Coordinate Reference System.
   • you can now use the one or the other at your will, as required and appropriate.