postgis 向量切片

文章主要來自於該篇文章，對部落格中的內容進行了注釋及部分修改，幫助理解。

#bbox函式，該函式要先在資料庫中執行才不會報錯，將文中的注釋去掉

create or replace function bbox(x integer, y integer, zoom integer)
returns geometry as
$body$
declare
max numeric := 6378137 * pi();#墨卡托投影的最大範圍
res numeric := max * 2 / 2^zoom;#墨卡托投影下乙個瓦片代表的實際距離
bbox geometry;
begin
根據傳入的行列號計算出該瓦片對應的墨卡托投影範圍並轉換到84座標系下，
返回的範圍是84座標系下的範圍--經緯度
return st_transform(st_makeenvelope(
-max + (x * res),
max - (y * res),
-max + (x * res) + res,
max - (y * res) - res,
3857),4326);
end;
$body$
language plpgsql immutable;

import psycopg2

#配置資料庫字典
dbconfig=
#該函式的作用就是根據表名，動態生成sql語句,在已知表結構的情況下，完全不需要這個函式。
def createsql(table_name,x,y,z):
geocolumn = "(select f_geometry_column from geometry_columns where f_table_name= '" + table_name + "' )"
#由於pg_attribute中沒有與表名直接對應的字段，先從pg_class中選出c表，得到其c.oid。再用a.attrelid判斷，即可得到c表
#實際上第一句sql做的事情就是將空間欄位列和屬性欄位列分開儲存f_geometry_column，和name列中
condition=" pg_class as c,pg_attribute as a  where c.relname = '" + table_name + "' and a.attrelid = c.oid and a.attnum>0"
sql1 = "select "+geocolumn+",replace(replace(string_agg(a.attname,','),"+geocolumn+"||',',''),"+geocolumn+"||'','') as name from" +condition
print(sql1)
conn = psycopg2.connect(database=dbconfig["database"], user=dbconfig["user"], password=dbconfig["password"], host=dbconfig["host"], port=dbconfig["port"])
cur = conn.cursor()
cur.execute(sql1)
keydata = cur.fetchall()
geom_column=keydata[0][0]
other_fields=keydata[0][1]
cur.close()
conn.close()
#select 前面就是正常的切片，from 後面的select在選擇空間屬性與普通屬性列，where條件則是先判斷是否相交，如果相交則返回交集
sql2="select st_asmvt(q, '"+table_name+"', 4096, 'geometry') from (select "+other_fields+"st_asmvtgeom("+geom_column+",bbox("+str(x)+", "+str(y)+" ,"+str(z)+") , 4096,256, true) as geometry  from "+table_name+" where "+geom_column+" && bbox("+str(x)+", "+str(y)+" ,"+str(z)+")  and st_intersects("+geom_column+", bbox("+str(x)+", "+str(y)+" ,"+str(z)+") )) q;"
return sql2;
#執行sql,獲得mvt的16進製制檔案
def getpbffile(table_name,x,y,z):
sql=createsql(table_name,x,y,z)
conn = psycopg2.connect(database=dbconfig["database"], user=dbconfig["user"], password=dbconfig["password"], host=dbconfig["host"], port=dbconfig["port"])
cur = conn.cursor()
cur.execute(sql)
print(sql)
keydata = cur.fetchall()
cur.close()
conn.close()
with open(str(z)+"_"+str(x)+"_"+str(y)+'.pbf', 'wb') as f:
f.write(keydata[0][0])
return keydata[0][0]
#x,y,z對應的是向量切片的行號，列號，級別
tile=getpbffile("cs",54605,26956,16)
print(tile)

postgis計算向量切片

沒寫錯，是使用postgis計算出來向量切片。在這之前先準備乙個資料 乙個gis資料表 本例中資料為一百萬的點資料，座標 4326 並在表中新增x,y欄位，方便後面的資料篩選。sql中用到了 st asmvt和st asmvtgeom。本文中建立向量切片很簡單,就是使用下方的乙個sql,執行結果如下...

gis 向量切片讀取 GIS向量切片演算法

structtilestructure public classvectortiletool listtiles publicvectortiletool public bool seprateshplayer string sourcepath,string resultfolder,intlev...

關於向量切片

憑回憶記一點 geoserver想切mapbox對應的向量切片的話，需要安裝對應的geoserver版本的外掛程式。ps 上傳不同格式的樣式檔案的時候也需要在這裡進行相應外掛程式的安裝 geoserver進行向量切片的時候，有好幾種可選的格式 geojson，topojson 記不清怎麼拼 以及對應...