项目需求，现有一批shp和tif数据，前端想要预览，但是没有服务发布环境；giser给他解决。
不废话、上代码。
1、搭建Python Web服务
我选择的是tornado框架，实现对外请求接口
import os
import sys
import tornado.web
import tornado.gen as gen
class MapnikTileServerHandler(tornado.web.RequestHandler):
@gen.coroutine
def get(self, rid, bsm, col, row, level):
self.set_header("Content-type", "image/png")
self.set_header('Access-Control-Allow-Origin', '*')
print('===========MapnikTileServerTileHandler============Starting...')
try:
rootPath = "{}/{}".format(cfg.TemDir, rid)
if not os.path.exists(rootPath):
self.write(u"无数据源，请检查。")
return None
mp = MapnikHandlerClass()
jsonFile_local = "{}/cof.json".format(rootPath)
with open(jsonFile_local, 'r') as fp:
key = json.loads(fp.read())
localPaths = json.loads(mp.b64decode(bsm))
pngFile = mp.main(rid, localPaths, key, float(col), float(row), level)
with open(pngFile, 'rb') as pic:
pics = pic.read()
self.write(pics)
except Exception as err:
log.writeLogWithTime(rid, str(err))
print (err)
finally:
print('==================================================End')2、引入mapnik
实现简单样式渲染和数据转换
def main(self, rid, localPaths, key,col=-1,row=-1,level=-1):
filePath = "{}/{}/L{}/{}_{}.png".format(cfg.MxdDir, rid, level, row, col)
if os.path.exists(filePath):
return filePath
datatype = key["dataType"].lower()
polygon_color = str(key["polygon_color"]) if "polygon_color" in key else "#F5F9F9"
line_color = str(key["line_color"]) if "line_color" in key else "#FF1400"
line_width = key["line_width"] if "line_width" in key else "0.5"
m = mapnik.Map(256, 256)
# m.background = mapnik.Color('steelblue')
# localPaths = ["D:\Data\shp/4407042014JBNTBHPK_ZW2_V2014.shp"]
box_main = None
for localPath in localPaths:
print (localPath)
name = localPath.replace('\\','/').split('/')[-1].split('.')[0]
print (name)
lyr = mapnik.Layer("{}_{}".format(name, datatype))
if datatype == "shp":
localPath = self.transform4326(localPath)
lyr.datasource = mapnik.Shapefile(file=localPath)
style_value = self.set_style_polygon(polygon_color, line_color, line_width)
style_name = "{}_Style".format(name)
m.append_style(style_name, style_value)
lyr.styles.append(style_name)
srs = lyr.srs
m.layers.append(lyr)
elif datatype == "tif":
localPath = self.transform4326(localPath)
localPath = self.create_map(localPath)
# lyr.datasource = mapnik.Gdal(file=localPath)
mapnik.load_map(m, localPath)
else:
return None
if "-1" in [str(row), str(col), str(level)]:
m.zoom_all()
savePath = "{}/{}/{}.png".format(cfg.TemDir, rid, uuid.uuid4())
mapnik.render_to_file(m, savePath)
else:
savePath = self.getTile(rid, m, col, row, level) #if datatype == "shp" else self.makeTile(rid, m, col, row, level,box_main)
print (savePath)
return savePath3、基于GDAL
实现数据转换，任何数据进来都按内部坐标系进行统一
def transform4326(self,input_lyr):
file_path, file_name = os.path.split(input_lyr)
out_lyr = "{}/copy_{}".format(file_path, file_name)
if os.path.exists(out_lyr):
return out_lyr
srs_proj4 = osr.SpatialReference()
srs_proj4.ImportFromEPSG(4326)
if input_lyr.endswith('.tif'):
# ds = gdal.Warp(out_lyr, input_lyr, srcSRS='EPSG:{}'.format(wkid), dstSRS=srs_proj4)
gdal.Warp(out_lyr, input_lyr, dstSRS=srs_proj4)
elif input_lyr.endswith('.shp'):
gdal.VectorTranslate(out_lyr, input_lyr, dstSRS=srs_proj4, reproject=True)
return out_lyr4、切片算法
根据前端传过来的行列号和层级进行切片计算，然后切成地图瓦片
def getTile(self,rid,m,col,row,level):
tiledirPath = "{}/{}".format(cfg.MxdDir, rid)
pathroot = "{}/L{}".format(tiledirPath, level)
imgname = "{}_{}.png".format(row, col)
filePath = "{}/{}".format(pathroot,imgname)
if os.path.exists(filePath):
return filePath
if not os.path.exists(pathroot):
os.makedirs(pathroot)
tilesize = 256
# m = mapnik.Map(tilesize, tilesize)
tileInfo = self.getResolution("4326")[str(level)]
step = float(tileInfo["resolution"]) * tilesize
origin = cfg.GetOrigin("4326")
minx = origin["x"] + step * float(col)
miny = origin["y"] - step * (float(row) + 1)
maxx = minx + step
maxy = miny + step
print ([minx,miny,maxx,maxy])
m.zoom_to_box(mapnik.Box2d(minx,miny,maxx,maxy))
mapnik.render_to_file(m, filePath)
return filePath5、整体实现
class MapnikHandlerClass():
def create_map(self,localPath):
file_path, file_name = os.path.split(localPath)
savePath = "{}/map_{}.xml".format(file_path, file_name.split('.')[0])
if os.path.exists(savePath):
return savePath
str_proj = "+proj=longlat +datum=WGS84 +no_defs +type=crs"
doc = Document()
map = doc.createElement('Map')
map.setAttribute('srs', str_proj)
doc.appendChild(map)
s = doc.createElement("Style")
s.setAttribute("name", "color relief style")
map.appendChild(s)
rule = doc.createElement("Rule")
s.appendChild(rule)
raster = doc.createElement("RasterSymbolizer")
raster.setAttribute("mode", "normal")
rule.appendChild(raster)
lyr = doc.createElement("Layer")
lyr.setAttribute("name", "color relief")
map.appendChild(lyr)
s_name = doc.createElement("StyleName")
text = doc.createTextNode("color relief style")
s_name.appendChild(text)
lyr.appendChild(s_name)
ds = doc.createElement("Datasource")
ty = doc.createElement("Parameter")
ty.setAttribute("name", "type")
p_ty = doc.createTextNode("gdal")
ty.appendChild(p_ty)
ds.appendChild(ty)
fe = doc.createElement("Parameter")
fe.setAttribute("name", "file")
p_fe = doc.createTextNode(localPath)
fe.appendChild(p_fe)
ds.appendChild(fe)
lyr.appendChild(ds)
with open(savePath, "w+") as f:
# 文件转码decode是为了保证输出的是UTF-8格式
f.write(doc.toprettyxml(encoding="UTF-8").decode("UTF-8"))
f.close()
return savePath
def transform4326(self,input_lyr):
file_path, file_name = os.path.split(input_lyr)
out_lyr = "{}/copy_{}".format(file_path, file_name)
if os.path.exists(out_lyr):
return out_lyr
srs_proj4 = osr.SpatialReference()
srs_proj4.ImportFromEPSG(4326)
if input_lyr.endswith('.tif'):
# ds = gdal.Warp(out_lyr, input_lyr, srcSRS='EPSG:{}'.format(wkid), dstSRS=srs_proj4)
gdal.Warp(out_lyr, input_lyr, dstSRS=srs_proj4)
elif input_lyr.endswith('.shp'):
gdal.VectorTranslate(out_lyr, input_lyr, dstSRS=srs_proj4, reproject=True)
return out_lyr
def b64decode(self,value):
u"""
base64解码
:param value: 待解码字符串
:return: 解码后的字符串
"""
return base64.b64decode(value + '=' * 3).decode('utf-8')
def b64encode(self,s):
u"""
base64编码
:param s: 待编码字符串
:return: 编码后的字符串
"""
return base64.b64encode(s.encode('utf-8')).decode('utf-8').replace('=', '')
def set_PolygonSymbol(self,rule, color='#F5F9F9'):
rule.symbols.append(mapnik.PolygonSymbolizer(mapnik.Color(color)))
def set_LineSymbol(self,rule, color='#FF1400',width=0.5):
rule.symbols.append(mapnik.LineSymbolizer(mapnik.Color(color), float(width)))
def set_style_polygon(self, polygon_color, line_color, line_width):
style = mapnik.Style()
rule = mapnik.Rule()
self.set_PolygonSymbol(rule,polygon_color)
self.set_LineSymbol(rule,line_color,line_width)
style.rules.append(rule)
return style
def set_RasterSymbol(self):
style = mapnik.Style()
rule = mapnik.Rule()
r = mapnik.RasterSymbolizer()
r.colorizer = mapnik.RasterColorizer()
for value, color in [
(0, "#7FFFD4"),
(1, "#76EEC6"),
(2, "#66CDAA"),
(3, "#458B74"),
(4, "#DAA520"),
(5, "#EEDD82"),
(6, "#7FFFD4"),
(7, "#76EEC6"),
(8, "#66CDAA"),
(9, "#458B74")
]:
r.colorizer.add_stop(value, mapnik.Color(color))
rule.symbols.append(r)
style.rules.append(rule)
return style
def getlyrinfo(self, dic_key,dataType):
extent = None
for i in dic_key:
localPath = self.transform4326(i)
if dataType == "shp":
dataset = ogr.Open(localPath) # 读取SHP数据
lyr = dataset.GetLayer(0) # 获取图层
extent = lyr.GetExtent()
xmin = extent[0]
xmax = extent[1]
ymin = extent[2]
ymax = extent[3]
extent = {'xmin': xmin, 'xmax': xmax, 'ymin': ymin, 'ymax': ymax}
elif dataType == "tif":
dataset = gdal.Open(localPath)
RasterYSize = dataset.RasterYSize
RasterXSize = dataset.RasterXSize
transform = dataset.GetGeoTransform()
x0 = transform[0]
y0 = transform[3]
cx = transform[1]
cy = transform[5]
xmin = x0
ymin = y0 + RasterYSize * cy
xmax = x0 + (RasterXSize * cx)
ymax = y0
extent = {'xmin': xmin, 'xmax': xmax, 'ymin': ymin, 'ymax': ymax}
return extent
def main(self, rid, localPaths, key,col=-1,row=-1,level=-1):
filePath = "{}/{}/L{}/{}_{}.png".format(cfg.MxdDir, rid, level, row, col)
if os.path.exists(filePath):
return filePath
datatype = key["dataType"].lower()
polygon_color = str(key["polygon_color"]) if "polygon_color" in key else "#F5F9F9"
line_color = str(key["line_color"]) if "line_color" in key else "#FF1400"
line_width = key["line_width"] if "line_width" in key else "0.5"
m = mapnik.Map(256, 256)
# m.background = mapnik.Color('steelblue')
# localPaths = ["D:\Data\shp/4407042014JBNTBHPK_ZW2_V2014.shp"]
box_main = None
for localPath in localPaths:
print (localPath)
name = localPath.replace('\\','/').split('/')[-1].split('.')[0]
print (name)
lyr = mapnik.Layer("{}_{}".format(name, datatype))
if datatype == "shp":
localPath = self.transform4326(localPath)
lyr.datasource = mapnik.Shapefile(file=localPath)
style_value = self.set_style_polygon(polygon_color, line_color, line_width)
style_name = "{}_Style".format(name)
m.append_style(style_name, style_value)
lyr.styles.append(style_name)
srs = lyr.srs
m.layers.append(lyr)
elif datatype == "tif":
localPath = self.transform4326(localPath)
localPath = self.create_map(localPath)
# lyr.datasource = mapnik.Gdal(file=localPath)
mapnik.load_map(m, localPath)
else:
return None
if "-1" in [str(row), str(col), str(level)]:
m.zoom_all()
savePath = "{}/{}/{}.png".format(cfg.TemDir, rid, uuid.uuid4())
mapnik.render_to_file(m, savePath)
else:
savePath = self.getTile(rid, m, col, row, level) #if datatype == "shp" else self.makeTile(rid, m, col, row, level,box_main)
print (savePath)
return savePath
def getResolution(self,wkid):
dic = {
"4326":{ str(i):{"resolution":1.40625/pow(2,i), "Scale":279541132.014359*2/pow(2,i)} for i in range(30)}
}
# for i in range(20):
# dic["4326"][str(i)] =
return dic[wkid] if wkid in dic else None
def getTile(self,rid,m,col,row,level):
tiledirPath = "{}/{}".format(cfg.MxdDir, rid)
pathroot = "{}/L{}".format(tiledirPath, level)
imgname = "{}_{}.png".format(row, col)
filePath = "{}/{}".format(pathroot,imgname)
if os.path.exists(filePath):
return filePath
if not os.path.exists(pathroot):
os.makedirs(pathroot)
tilesize = 256
# m = mapnik.Map(tilesize, tilesize)
tileInfo = self.getResolution("4326")[str(level)]
step = float(tileInfo["resolution"]) * tilesize
origin = cfg.GetOrigin("4326")
minx = origin["x"] + step * float(col)
miny = origin["y"] - step * (float(row) + 1)
maxx = minx + step
maxy = miny + step
print ([minx,miny,maxx,maxy])
# if str(int(level)) == "0":
# m.zoom_all()
# else:
m.zoom_to_box(mapnik.Box2d(minx,miny,maxx,maxy))
mapnik.render_to_file(m, filePath)
return filePath
6、切片效果
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
