地理信息系统设计与开发（空间查询和创建缓冲区）.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 地理信息系统设计与开发（空间查询和创建缓冲区）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import ogr\n",
    "import osr\n",
    "import json\n",
    "import requests\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 设置输入、输出文件和缓冲区距离（与空间参考一致）"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputfn = 'test3.shp'\n",
    "outputBufferfn = 'testBuffer4.shp'\n",
    "bufferDist = 5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 获取输入数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "inputds = ogr.Open(inputfn)\n",
    "inputlyr = inputds.GetLayer()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 由json创建一个包含全中国区域的multipolygon，用于空间查询"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def createMultiPolygonFromMultiPolygonJson():\n",
    "    myMultipolygon = ogr.Geometry(ogr.wkbMultiPolygon)\n",
    "    ds=ogr.Open('china.json')\n",
    "    lyr=ds.GetLayer() \n",
    "    for feature in lyr:\n",
    "        geom=feature.geometry()\n",
    "        if(geom.GetGeometryName() == 'MULTIPOLYGON'):\n",
    "            for i in range(geom.GetGeometryCount()):\n",
    "                ring=geom.GetGeometryRef(i)\n",
    "                #print('ring:',ring)\n",
    "                myMultipolygon.AddGeometry(ring)\n",
    "        else:\n",
    "            myMultipolygon.AddGeometry(geom)\n",
    "    return myMultipolygon "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 进行空间查询并返回查询结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spatialFilter(inputlyr):\n",
    "    myMultiPolygon=createMultiPolygonFromMultiPolygonJson()\n",
    "    print(myMultiPolygon.GetGeometryCount())\n",
    "    print(\"before filter:\",inputlyr.GetFeatureCount())\n",
    "    inputlyr.SetSpatialFilter(myMultiPolygon)\n",
    "    print(\"after filter:\",inputlyr.GetFeatureCount())\n",
    "\n",
    "    return inputlyr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "288\n",
      "before filter: 14011\n",
      "after filter: 16\n"
     ]
    }
   ],
   "source": [
    "result=spatialFilter(inputlyr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 显示查询结果"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "127km WSW of Aksu, China\n",
      "127km WSW of Aksu, China\n",
      "66km ENE of Aksu, China\n",
      "59km E of Sary-Tash, Kyrgyzstan\n",
      "154km ENE of Luring, China\n",
      "61km NNW of Pubu, China\n",
      "87km ESE of Arzak, China\n",
      "31km E of Junlian, China\n",
      "6km SW of Changning, China\n",
      "99km E of Arzak, China\n",
      "22km NNW of Biruxong, China\n",
      "54km N of Ziro, India\n",
      "111km SSE of Kuqa, China\n",
      "3km NW of Fuji, China\n",
      "84km SSE of Kuqa, China\n",
      "229km SE of Kuqa, China\n"
     ]
    }
   ],
   "source": [
    "resultFeat = result.GetNextFeature()\n",
    "while resultFeat :\n",
    "   print(resultFeat.GetField('place')) \n",
    "   resultFeat = result.GetNextFeature()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 创建输出文件"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "shpdriver = ogr.GetDriverByName('ESRI Shapefile')\n",
    "if os.path.exists(outputBufferfn):\n",
    "    shpdriver.DeleteDataSource(outputBufferfn)\n",
    "outputBufferds = shpdriver.CreateDataSource(outputBufferfn)\n",
    "targetSR = osr.SpatialReference()\n",
    "targetSR.ImportFromEPSG(4326) #Geo WGS84\n",
    "bufferlyr = outputBufferds.CreateLayer(outputBufferfn,targetSR,geom_type=ogr.wkbPolygon)\n",
    "featureDefn = bufferlyr.GetLayerDefn()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 对查询结果建立缓冲区并输出为文件"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "for feature in inputlyr:\n",
    "    ingeom = feature.GetGeometryRef()\n",
    "    geomBuffer = ingeom.Buffer(bufferDist)\n",
    "\n",
    "    outFeature = ogr.Feature(featureDefn)\n",
    "    outFeature.SetGeometry(geomBuffer)\n",
    "    bufferlyr.CreateFeature(outFeature)\n",
    "    outFeature = None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "outputBufferds=None"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}