QGIS3.6在MAC中打开文件夹速度奇慢

忘记从QGIS3哪个版本开始,在MAC中打开目录对话框就速度奇慢,导致我一开始都是直接将文件的绝对路径直接键入进来,避免打开文件对话框。搜索发现,此bug已有用户汇报了,且给出了一个解决方法:

在系统偏好设置->安全性与隐私->隐私->辅助功能中,将QGIS3.6的选项给去掉即可。

参考:https://gis.stackexchange.com/questions/310280/qgis-very-slow-when-browsing-files-on-macos

升级greenplum 5.x到6.0beta3

之前的老集群,使用的是ubuntu 16.04下的greenplum 5.8,后来随系统升级到5.19版本。最近发现运行速度变慢(慢得过分,可以一个sql喝几杯茶的程度),而且经过和在另个一个新集群上的相同数据对比测试,发现速度差别至少有2个量级。

首先尝试将数据库备份,删除再恢复,发现速度改善很小。有了上个集群的安装经验,于是决定将操作系统升级,并且将greenplum也升级6.0版本(其master版本是7.0,但稳定性还是有问题,发现有几个sql查询异常)。

在升级过程中,遇到一个小坑:机器名不要命名为host01-1的形式,会导致greenplum将其判别为一类特殊的host。这个问题还是通过浏览gpinitsystem的代码才找出来的。

另外也再记录一下postgis的安装过程。由于libxerces-c-dev和gp-xerces的冲突,postgis安装时需要将raster功能禁止,在greenplum安装完成后执行:

source /usr/local/gdbp/greenplum_path.sh
git clone https://github.com/greenplum-db/geosptial.git
cd geospatial/postgis/build/postgis-2.1.5
./configure --with-pgconfig=$GPHOME/bin/pg_config --without-raster --without-topology --prefix=$GPHOME
make USE_PGXS=1 clean all
sudo make USE_PGXS=1 install 

这样在所有节点安装完成后即可开始greenplum的初始化脚本。

全国县市边界数据的处理(1M2017版)

全国地理信息资源目录服务系统上可以下载最新版本(2017年,1:1百万比例尺)的全国各要素基础数据,具体包括:

  • 共77幅DLG数据集,包括全国陆地范围以及台湾岛、海南岛、钓鱼岛、南海诸岛在内的主要岛屿及其临近海域。
  • 包括水系、居民地及设施、公路铁路、行政边界、地名等5个数据集12个数据层。

我们现在就需要其中的行政边界数据,本来以为下载后可以很方便的使用(有点想当然了),经过一番艰苦的处理后,数据终于是可用状态了。

在mac系统下进行了数据处理,使用了gdal/ogr、QGIS3,主要处理过程记录如下。

  1. 将下载下来的77副gdb数据,全部解压到一个目录中,比如1M2017目录,并在其中建立一个shp子目录。
  2. 在shell下,使用ogr2ogr命令,将所有gdb文件转换为一个大的shape文件,注意需要将其中的NAME字段编码转换为UTF8编码(官方是GBK编码)。
  3. 使用ogr2ogr,对合并后的BOUA(面状边界数据)文件根据PAC字段进行再次图形合并(因为默认的分幅会保留大量的正方形),生成output1.shp文件。
  4. output1.shp文件包含了大量了国外数据、海洋数据,这些数据可以通过ogr2ogr再次剔除,即设置PAC字段>100000且不等于250100,然后生成output2.shp文件。
  5. 仔细观察output2.shp文件,发现其中存在大量的拓扑错误(分幅边线的交叉遗留),在QGIS3中使用v.clean工具进行错误检查和处理,生成output3.shp文件。
  6. 对output3.shp文件进行人工错误检查(可借助QGIS3的图形检查工具check validity辅助检查),使用编辑工具(以及高级数字化工具)进行错误修复,主要使用节点编辑工具(删除节点)、多边形合并工具、多边形part删除工具等,最终生成全国(含港澳台)的县级边界数据集boua_cnty.shp。
  7. 使用QGIS3在boua_cnty数据中添加一个虚拟字段pac_city,并将其值定义为floor(PAC/100),然后对其执行dissolve操作,若需要将直辖市也当做一个独立城市边界,则需要单独处理一下四个直辖市的属性(将县区进行合并,主要是天津和重庆),形成全国市级边界数据集boua_city.shp文件。
for i in `ls -d *.gdb`;do ogr2ogr -f "ESRI Shapefile" -lco ENCODING=UTF-8  -append  shp $i;done
ogr2ogr -f "ESRI Shapefile" output1.shp -lco ENCODING=UTF-8 -dialect sqlite -sql "select st_union(geometry) as geometry,PAC,NAME from 'BOUA' group by PAC,NAME order by PAC,NAME" shp/BOUA.shp
ogr2ogr -f "ESRI Shapefile" output2.shp -lco ENCODING=UTF-8 -where "PAC>100000 and PAC<>250100" output1.shp

但这样生成的城市边界仅包括陆地边界,对于内陆型城市,在使用中问题不大,但对于沿海城市,则很多人类活动轨迹可能会超出陆地边界(比如舟山),因此还需要一个融合近海边界的城市边界。手头上有一个2012年版本的含海洋边界的城市边界数据,需要将其海洋边界数据融合到新的全国边界数据集中。

处理过程为:

  1. 对含有海洋边界的数据集bou2012,和最新全国城市边界数据集boua_city,执行difference操作,在QGIS3进行编辑,删除内陆边界中有差异的部分,在沿海边界删除海洋中的小岛(内多边形),同时添加一个虚拟字段pac_city(和boua_city中的字段相同),最终生成海洋城市边界数据集。
  2. 将boua_city数据集和上一步中生成的海洋城市边界数据集执行union操作。
  3. 对union后的数据集,根据pac_city字段执行dissolve操作。
  4. 人工检查最终生成的数据集。

因数据问题,此数据集中不含港澳台的海洋边界。

PostGIS:ST_xxx函数中ST有什么含义

此前先入为主,一直以为ST的含义是Spatial and Temporal。在课堂上,有学生问我这个含义,我以这个含义进行了解释。但心里感觉没有太大的把握,因此课后回来专门搜索了一下,发现网络上还有人问了这个问题:

https://stackoverflow.com/questions/7234679/what-is-st-in-postgis

从响应来看,最早ST的确是Spatial和Temporal的含义,但在PostGIS的官方文档中是作为Spatial Type进行解释的:

the standard spatial type (ST) prefix.

Window环境下qgis standalone程序开发环境设定

本文假定:
1. 采用qgis官方发布的安装程序(如QGIS-OSGeo4W-2.18.13-1-Setup-x86.exe安装包)
2. 使用默认的安装位置,即C:\Program Files\QGIS 2.18,如指定到其他目录,请修改对应的设定。
3. 假定示例数据目录为D:\QGIS\数据\

第一:设定环境变量

SET OSGEO4W_ROOT=C:\Program Files\QGIS 2.18
SET QGIS_PREFIX=%OSGEO4W_ROOT%\apps\qgis
SET PATH=%QGIS_PREFIX%\bin;%OSGEO4W_ROOT%\bin
SET PYTHONPATH=%QGIS_PREFIX%\python;%OSGEO4W_ROOT%\apps\Python27
SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python27
SET GDAL_DATA=%OSGEO4W_ROOT%\share\gdal

即在cmd下执行上述命令,若路径不一致,请自行修改。

第二:运行一个样例

# -*- coding: utf-8 -*-
import os,os.path,sys

from qgis.gui import *
from qgis.core import *
from PyQt4 import QtGui

app = QtGui.QApplication([])
QgsApplication.setPrefixPath(os.environ['QGIS_PREFIX'],True)
QgsApplication.initQgis()

main_win = QtGui.QMainWindow()
frame = QtGui.QFrame(main_win)
main_win.setCentralWidget(frame)
grid_layout = QtGui.QGridLayout(frame)

map_canvas=QgsMapCanvas()
grid_layout.addWidget(map_canvas)
map_canvas.setCanvasColor(QtGui.QColor(255,255,255))
layer=QgsVectorLayer(ur'D:\qgis\数据\bou2_4p.shp','border','ogr')

QgsMapLayerRegistry.instance().addMapLayer(layer)
canvas_layer=QgsMapCanvasLayer(layer)
map_canvas.setLayerSet([canvas_layer])
map_canvas.zoomToFullExtent()

main_win.show()
# Need the following statement if running as a script
app.exec_()

可在前述cmd环境下运行python,将此代码粘贴进来运行。

第三,说明:
layer=QgsVectorLayer(ur'D:\qgis\数据\bou2_4p.shp','border','ogr')
这个语句,如目录名称中不含中文,则路径前的u可以移除。当然,也可以这样写:
layer=QgsVectorLayer(u'D:\\qgis\\数据\\bou2_4p.shp','border','ogr')
同时,建议在python下逐行执行来理解代码。

QGIS中TimeManager的一个问题修复

在QGIS中使用TimeManager时遇到一个问题,即满足时空要素要求后,时空数据展示不出来。点击数据的layer查看发现其对应的sql语句大概是:
dbname='citybigdata' host=localhost port=5432 sslmode=disable key='id' srid=4326 type=Point table="qqheat"."milocal" (geom) sql=cast("ts" as character) < '2016-11-09 21:05:00' AND cast("ts" as character) >= '2016-11-09 21:00:00'
将类似的语句直接在psql中发现执行的时候有问题,无法返回结果。但将语句修改一下即可解决:
dbname='citybigdata' host=localhost port=5432 sslmode=disable key='id' srid=4326 type=Point table="qqheat"."milocal" (geom) sql=cast("ts" as text) < '2016-11-09 21:05:00' AND cast("ts" as text) >= '2016-11-09 21:00:00'
即将其中的character改成text即可。
为了使插件工作正常,就需要修改timemanager的代码,使用grep工具找到该代码的位置:

cd .qgis2/python/plugins/timemanger
grep -rnw ‘.’ -e ‘character’

发现是如下代码的问题:

./.git/hooks/pre-commit.sample:26:# printable range starts at the space character and ends with tilde.
./CONTRIBUTING.md:13:* Go to http://www.loc.gov/standards/iso639-2/php/code_list.php and find the 2-character ISO 639-1
./docs/Doxyfile:586:# This tag can be used to specify the character encoding of the source files
./query_builder.py:9:STRINGCAST_FORMAT = ‘cast(“{}” as character) {} \'{}\’ AND cast(“{}” as character) >= \'{}\’ ‘
./query_builder.py:56: return ‘ cast(“{}” as character) LIKE \’%BC\”.format(attr)
./query_builder.py:63: return ‘ cast(“{}” as character) LIKE \’%AD\”.format(attr)
./query_builder.py:79: return ” ‘{}’ {} cast(\”{}\” as character) “.format(val, comparison, col)
./query_builder.py:87: return ” ‘{}’ {} cast(\”{}\” as character) “.format(val, comparison, col)

将具体的代码修改,然后重启qgis,即可工作正常。

使用 nginx+Tilestache 构建瓦片地图

之前有介绍 “使用tilestache构建瓦片地图服务器”,在我们的服务器上使用一段时间后,发现用其作为server,稳定性有一定问题。因此就考虑使用nginx直接提供静态地图文件,而只使用tilestache来生产瓦片。tilestache的瓦片缓存是按照一定的目录规则来实现的,即在url中表现为/layername/level/x/y.png,而在磁盘上存储的方式为/layername/level/umask/x/umask/y.png,所以需要进行rewrite转向。网上搜索后发现已经有人解决了这个问题:

http://andyfiedler.com/blog/using-nginx-to-proxy-tilestache-180/

 

在服务器上,按照上述的修改,即可满足实际需求。

 

使用tilestache构建瓦片地图服务器

Tilestache是一款开源的瓦片地图服务器,其安装和使用较简单,但官方的文档感觉还是偏简单了些。在这儿记录一下最近安装和使用tilestache的方法。

1、 安装

在Ubuntu服务器上安装tilestache还是比较简单的:

sudo apt-get install tilestache python-pil python-gdal

tilestache默认并没有绑定python-pil和python-gdal两个库,若开始只是安装了tilestache,后面可能可能会出现返回空白地图的问题。

2. 使用

tilestache使用时要构建一个cfg的配置文件,同时若采用mapnik来进行渲染,还必须提供每个对应图层的渲染配置文件。

一个简要的配置文件示例如下:

{
  "cache": {
    "name": "Disk",
    "path": "cache/"
  },
  "layers":
  {
   "landuse2005":
        {
        "provider":{"name":"mapnik","mapfile":"landuse2005.xml"}, "projection": "spherical mercator",
        "preview":   {        "lat": 33.4,       "lon": 97.3,        "zoom": 15      }
        },
    "roads":
   {
        "provider":{"name":"proxy","url":"http://tile.openstreetmap.org/{Z}/{X}/{Y}.png"},
	"preview":   {        "lat": 33.4,       "lon": 97.3,        "zoom": 15      }
    }
  }
}

示例中有两个图层,landuse2005采用mapnik进行渲染,roads采用OSM的瓦片并直接进行代理。mapnik的渲染文件入门不易,可以采用tilemill来生成处理。

3. 预制瓦片

上述图层的landuse2005应该预先生成,用户访问时才能提高速度。tilestache同时也提供了一个程序来生成:

tilestache-seed -b 31.8 89.3 36.3 102.3 -c tilestache.cfg -l landuse2005 10 11 12 13 14 15

4. 启动

最简单的方法就是使用tilestache-server来启动对应的服务:

tilestache-server -c tilestache.cfg -i 192.168.13.13

参数比较简单,-i就是要绑定的网络接口,然后用浏览器访问:

http://192.168.13.13:8080/landuse2005/preview.html

 

 

利用shell命令合并分县土地利用数据为一个shape File

为了制图服务的需求,需要将分县的土地利用arcinfo coverage数据合并为一个shape File。
用ogrinfo查看一个县的数据简要信息:
ogrinfo -so id04632323 PAL
Had to open data source read-only.
INFO: Open of `id04632323'
using driver `AVCBin' successful.

Layer name: PAL
Geometry: Polygon
Feature Count: 4065
Extent: (-831725.500000, 2525475.000000) - (-704241.375000, 2614980.250000)
Layer SRS WKT:
PROJCS["unnamed",
GEOGCS["Unknown datum based upon the Krassowsky 1940 ellipsoid",
DATUM["Not_specified_based_on_Krassowsky_1940_ellipsoid",
SPHEROID["Krassowsky 1940",6378245,298.3,
AUTHORITY["EPSG","7024"]],
AUTHORITY["EPSG","6024"]],
PRIMEM["Greenwich",0,
AUTHORITY["EPSG","8901"]],
UNIT["degree",0.0174532925199433,
AUTHORITY["EPSG","9122"]],
AUTHORITY["EPSG","4024"]],
PROJECTION["Albers_Conic_Equal_Area"],
PARAMETER["standard_parallel_1",25],
PARAMETER["standard_parallel_2",47],
PARAMETER["latitude_of_center",12],
PARAMETER["longitude_of_center",110],
PARAMETER["false_easting",0],
PARAMETER["false_northing",0],
UNIT["METERS",1]]
ArcIds: IntegerList (0.0)
AREA: Real (12.3)
PERIMETER: Real (12.3)
ID04632323#: Integer (5.0)
ID04632323-ID: Integer (5.0)

用QGIS打开后,可以发现其ID04632323-ID字段是土地利用的分类代码。而要将这些分县的数据合并,则需要利用ogr2ogr来合并生成。
主要问题:
1. 需要重命名字段,如将这个分类代码重命名为一个共用的名称,比如code。
2. 目录名和分类代码字段名相关,但需要进行大小写转换。
3. 还需要将每个县的目录名称保存到对应的数据库中。

在bash中,大小写转换可以用这个函数,将变量$i转换为大写:
${i^^}

ogr2ogr可以用sql选项来提取不同的字段并重命名:
select AREA as area,PERIMETER AS perimeter,ID04632323-ID as code from PAL
但在测试中,发现其不认这种写法,需要将其用单引号来括起来:
select AREA as area,PERIMETER AS perimeter,‘ID04632323-ID’ as code from PAL

而目录名称也可以通过变量的方式传递到sql中来。

最终命令脚本,要注意其中的引号和转义:
for i in id*;do
ogr2ogr -update -append -sql "select AREA as area,PERIMETER AS perimeter,'${i^^}-ID' as code,\"$i\" as county FROM PAL" landuse.shp $i;
done