利用python加gdal将shp文件导入postgis数据库相关问题

---

前言

本人是某高校在读学术垃圾，以下经验仅供参考，如有网络高人不吝赐教，小弟感激不尽。最近本人开始肝毕设，导师要求建数据库以便存放毕设中原始数据及成果数据（感觉是老师怕我论文内容不够），于是我便想起了gdal。查阅资料后发现，gdal支持比较好的数据库有postgis和Oracle的空间数据库，对比后个人感觉还是postgis简单一些，于是便试着学一学。postgis安装教程内网很容易找到，只是安装过程中，因人而异，会有些稀奇古怪的bug。一番折腾以后，我总算是部署好数据库了，但麻烦事才刚刚开始。泪目，加油我能毕业的！！！

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、如何利用gdal将shp文件导入到postgis数据库？

官方文档给的用法很简单，用以下命令即可：

// ogr2ogr表示要使用的程序，overwrite表示覆盖写入
//f我也忘了是什么
//PG后面紧跟的字符串是你的数据的相关访问信息
//host是服务器可以填地址（没试过）或直接localhost（表示本地）
//dbname为你要导入的数据库名，其他都很好懂
//最后是文件路径
//更多相关参数具体见官方文档
ogr2ogr -overwrite -f "PostgreSQL" PG:"host=localhost dbname=mydb user=postgres password=xxxxx" F:\xxx\test\quj.shp

但是这样在python中是运行不了的。于是，我又在网上找方法，把某度的肠子都翻出来了，只找到一篇文章，貌似就是我想要的：
链接: https://zhuanlan.zhihu.com/p/160166790
文章告诉我的这样用：

// An highlighted block
import os
os.system('ogr2ogr '+'-overwrite '+'-f '+'"'+"PostgreSQL"+'"'+' PG:'+'"'+"host=localhost user=postgres dbname=mydb password=xxxx"+'"'+' '+'"'+"F:\xxx\test\quj.shp"+'"')

但是还是不行，而且运行完只会返回数字1。我：啥意思？
于是我提问了文章的作者，作者建议让我现在CMD上试试。（这里感谢知乎用户：一个有趣的灵魂W，热情帮助。）
我这才意识到，原来gdal官方给的是命令行命令啊？于是跑到CMD运行了最开始的那段代码，结果报错：
‘ogr2ogr’ 不是内部或外部命令，也不是可运行的程序或批处理文件。
后来才发现，ogr2ogr是可执行程序。所以最后把ogr2ogr的路径也带上就能成功运行。
python中的正确调用：

order1='D:\\anaconda3\\Lib\\site-packages\\osgeo\\ogr2ogr '+'-overwrite '+'-f '+'"'+"PostgreSQL"+'"'+' PG:'+'"'+"host=localhost user=postgres dbname=postgis_31_sample password=xxxx"+'"'+' '+'"'+"F:\\qinshigou\\test\\szt1.shp"+'"'
os.system(order1)
//网上说可以用os.system('cd /路径 && 后续命令')的方法运行，但不知道为啥我这里cd不了。所以直接把两段命令拼在一起了。

二、结果

1.结果

可以看到导入是成功的，但还有别的报错。

2.报错

如下问题：

ERROR 1: PROJ: proj_identify: D:\postgresql\10\share\contrib\postgis-3.1\proj\proj.db lacks DATABASE.LAYOUT.VERSION.MAJOR / DATABASE.LAYOUT.VERSION.MINOR metadata. It comes from another PROJ installation.

意思是说，postgis的投影库中找不到对应的投影。

---

对上述报错一开始我没太在意，一是在数据库中看数据没问题，二是看不大懂在说什么玩意。凭借我的渣滓英语勉强知道在proj.db里匹配不到啥东西，起初还以为是我安装的时候又做错了什么。

后来试了导入再导出2份数据，才明白这个报错是啥意思。
第一份数据是wgs84坐标，发现导出的跟导入的坐标系变了：
从GCS_WGS_1984变成了WGS 84，对比后发现，只是名字变了，这可以忽略。
第二份数据是wgs84下的投影坐标系，从库里导出后坐标直接没了：
原数据是Albers_Conical_Equal_Area，阿尔伯斯等面积圆锥投影。

经过检测proj.db发现，在Admin4中应该对应的是spatial_ref_sys这张表，postgis里貌似没有这个投影。
于是把第二份数据换成了UTM投影后导入再导出就正常了，同时上述报错也没了。

总结

gdal官方给的postgis导入shp的命令实际上是命令行命令，需要利用ogr2ogr这个可执行程序来做。在Python中需要调用os.system(字符串)来执行该命令行命令。

需要注意的是：
1.要cd到ogr2ogr的目录下执行命令。
2.导入数据前先看看postgis支持哪些坐标系，尽量采用知名的，国外的WGS84/UTM，国内的北京54，西安80和大地2000都行。不要像我一样用一些比较奇怪的坐标系，以免数据出问题。

gdal还支持栅格数据导出postgis，但不支持导入（没记错的话），下一步我试试栅格数据的导出怎么用Python+gdal来做。
官方画饼说未来导入导出都支持，希望他们早点做出来吧。

后续测试发现，其实坐标系的那个报错即使换成了UTM或WGS84、北京54等还是会出现，只是坐标系不再丢失。猜想应该是Arcgis的坐标名称跟postgis 的不同吧。

勘误（2021-08-11）

ERROR 1: PROJ: proj_identify: D:\postgresql\10\share\contrib\postgis-3.1\proj\proj.db lacks DATABASE.LAYOUT.VERSION.MAJOR / DATABASE.LAYOUT.VERSION.MINOR metadata. It comes from another PROJ installation.

上述报错是我当时误解了，“It comes from another PROJ installation.”直译为“这来自于其他PROJ安装”。翻译成人话就是调用的proj.db跟gdal不匹配。
因为近来使用gdal处理数据时经常出现该报错，严重的情况还导致新建数据时坐标导不进去。多番调试没解决，最后受到stackexchange上的帖子的启发（这网站居然不需要魔法上网，还真是神奇。），最终我明白了问题出在哪。
其中这个帖子：（https://gis.stackexchange.com/questions/378463/cannot-find-proj-db-and-error-1-proj-proj-create-from-database-cannot-buil）中有人说：他自己遇到了帖子中同样的报错（“Cannot find proj.db” and “ERROR 1: PROJ: proj_create_from_database: cannot build geodeticCRS 4326” errors），同时他安装了postgis。postgis在系统中创建了proj.db的环境变量，使得gdal引用了postgis的proj.db，而非gdal自己的，这样就会导致出错。

看完上述发言，我瞬间懂了，原因是：gdal和postgis给proj.db设置的环境变量的名字是一样的。我打开环境变量表看，发现只有postgis的（估计是我先装的postgis，后来的gdal加不进去，因为重名）。然后我运行gdal时，gdal这憨憨就以为postgis的proj.db是它的。结果当然就是会报错。

意识到这个问题后，我在代码开头

import os
os.environ['PROJ_LIB'] = 'C：user\proj'\\假设gdal对应的proj.db在这个文件夹下

现在再运行gdal就不会报上述错误，数据的坐标系也不会出些莫名其妙的问题了。