postgis

在gis中前端通常采用cesuim加载terrian文件的形式来计算海拔（高程），但是避免不了数据量大的情况会出现页面卡顿，于是想办法在后端尝试计算。网上查了一些资料也是比较零散，于是具有缝合怪属性的我开始缝缝补补，我们不生产代码，只做代码搬运工。开始采用的是python的gdal库进行计算，可以实现而且响应也还可以，6000条的数据在3s左右，不过身为程序员的“苦逼们”不是在优化就是在优化的路上，这里也是查到了一些postgre中的一个扩展postgis，它可以将栅格（.tif文件）数据转成sql并通过函数直接计算，尝试了一下效果显著，这不得让我在简历上浓墨重彩一下。因为是第一次接触，里面的函数还不是很清楚，都是大佬们的智慧结晶。下面是具体的搭建方式以及测试用例

安装postGIS

#创建数据卷
mkdir -p /docker/postgis/data
#拉取镜像
docker pull postgis/postgis
#运行容器
docker run --name postgres --restart always \
-e POSTGRES_PASSWORD=123456 \
-v /docker/postgis/data:/var/lib/postgresql/data \
--privileged=true \
-p 5432:5432 \
-d postgres:13

修改密码
docker exec -it  ps命令查出来的id bash
su postgres
psql -U postgres
Alter user postgres with password '123456';

创建拓展
CREATE EXTENSION postgis;
CREATE EXTENSION pgrouting;
CREATE EXTENSION postgis_topology;
CREATE EXTENSION fuzzystrmatch;
CREATE EXTENSION postgis_tiger_geocoder;
CREATE EXTENSION address_standardizert;

安装raster2pgsql
apt-get update
apt-get install postgis

将tif转sql
raster2pgsql -s 4326 -I -C -M -F -t 256x256 *.tif public.jsdemelevation > elev.sql
导入数据
psql -d postgres -f elev.sql -U postgres -W

-- 查询一个点高程
SELECT 
	ST_Value(rast, ST_SetSRID(ST_MakePoint(118.9620236, 32.23440476), 4326)) val
FROM
	jsdemelevation
WHERE
	ST_Intersects(rast, ST_SetSRID(ST_MakePoint(118.9620236, 32.23440476), 4326))



-- 查询指定经纬度的高程数据（假设查询的经纬度为lng_value和lat_value）
SELECT
    ST_Value(rast, geom) AS val
FROM
    jsdemelevation,
    (VALUES
        (ST_SetSRID(ST_MakePoint(118.9620236, 32.23440476), 4326)),
        (ST_SetSRID(ST_MakePoint(lon2, lat2), 4326)), -- 第二个点的经纬度
        (ST_SetSRID(ST_MakePoint(lon3, lat3), 4326))  -- 第三个点的经纬度
    ) AS points(geom)
WHERE
    ST_Intersects(rast, geom);
    
-- 查询指定空间范围内某条线上各个点位置对应的高程值
SELECT 
	ST_Value(dm.rast, ps.geom) val
FROM
	jsdemelevation AS dm, (
		SELECT 
			(ST_DumpPoints(geom)).geom
		FROM
			a_line_data
		WHERE 
			ST_Intersects(geom, ST_MakeEnvelope(
			118.9620236, 32.23440476, 4326))
	) AS ps
WHERE
	ST_Intersects(dm.rast, ps.geom)

测试用例，10000条数据响应稳定在1s内，不建议用java，java mybatis中标签进行拼接会耗时大量时间，可以采用python字符串拼接的方式

#!/usr/bin/env python3
# coding: utf-8
from flask import Flask, request, jsonify, Response
import psycopg2
import json

app = Flask(__name__)
ALLOWED_EXTENSIONS = {'json'}

# 数据库连接参数
dbname = 'postgres'
user = 'postgres'
password = 'password'
host = 'ip'
port = '5432'

# 连接到数据库
conn = psycopg2.connect(dbname=dbname, user=user, password=password, host=host, port=port)
cur = conn.cursor()


def allowed_file(filename):
    return '.' in filename and filename.rsplit('.', 1)[1].lower() in ALLOWED_EXTENSIONS


@app.route('/elevation', methods=['POST'])
def get_elevation():
    if 'file' not in request.files:
        return Response({'error': 'Invalid request'}, status=400, mimetype='application/json')
    file = request.files['file']
    if file.filename == '':
        return Response({'error': 'Invalid request'}, status=400, mimetype='application/json')
    if file and allowed_file(file.filename):
        try:
            # 构建临时表的值部分
            values_part = ''
            json_stream = file.stream
            data_list = json.load(json_stream)
            for idx, item in enumerate(data_list):
                # 提取经纬度信息
                latitude = item['latitude']
                longitude = item['longitude']
                # 构建每个点的值部分，并进行拼接
                if idx == 0:
                    values_part += f"(ST_SetSRID(ST_MakePoint({longitude}, {latitude}), 4326))"
                else:
                    values_part += f",\n(ST_SetSRID(ST_MakePoint({longitude}, {latitude}), 4326))"
            # SQL查询模板
            sql_query = f"""
            SELECT
                ST_Value(rast, geom) AS val
            FROM
                jsdemelevation,
                (VALUES
                    {values_part}
                ) AS points(geom)
            WHERE
                ST_Intersects(rast, geom);"""
            # 执行查询
            cur.execute(sql_query)

            # 获取所有结果
            results = cur.fetchall()
            return jsonify({'code': 200, 'data': results, 'msg': 'success'})
        except Exception as e:
            print(e)
            # 如果发生任何异常，返回状态码 400
            return Response({'error': 'Invalid request'}, status=400, mimetype='application/json')
    else:
        return jsonify({'code': 400, 'data': None, 'msg': 'No file part'})


if __name__ == '__main__':
    app.run(debug=True, port=8000)