历史文章附件列表/（数据科学学习手札75）基于geopandas的空间数据分析——坐标参考系篇/（数据科学学习手札75）基于geopandas的空间数据分析——坐标参考系篇.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore') # 屏蔽所有警告信息"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3 geopandas中的坐标参考系管理"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Projected CRS: +proj=tmerc +lat_0=0 +lon_0=108 +k=1 +x_0=500000 + ...>\n",
       "Name: unknown\n",
       "Axis Info [cartesian]:\n",
       "- E[east]: Easting (metre)\n",
       "- N[north]: Northing (metre)\n",
       "Area of Use:\n",
       "- undefined\n",
       "Coordinate Operation:\n",
       "- name: unknown\n",
       "- method: Transverse Mercator\n",
       "Datum: Unknown based on IAU 1976 ellipsoid\n",
       "- Ellipsoid: IAU 1976\n",
       "- Prime Meridian: Greenwich"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import pyproj\n",
    "\n",
    "pyproj.CRS.from_user_input('+proj=tmerc +lat_0=0 +lon_0=108 +k=1 +x_0=500000 +y_0=0 +ellps=IAU76 +units=m +no_defs')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Projected CRS: EPSG:2381>\n",
       "Name: Xian 1980 / 3-degree Gauss-Kruger CM 108E\n",
       "Axis Info [cartesian]:\n",
       "- X[north]: Northing (metre)\n",
       "- Y[east]: Easting (metre)\n",
       "Area of Use:\n",
       "- name: China - 106.5°E to 109.5°E onshore\n",
       "- bounds: (106.5, 18.19, 109.5, 42.47)\n",
       "Coordinate Operation:\n",
       "- name: 3-degree Gauss-Kruger CM 108E\n",
       "- method: Transverse Mercator\n",
       "Datum: Xian 1980\n",
       "- Ellipsoid: IAG 1975\n",
       "- Prime Meridian: Greenwich"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyproj.CRS.from_user_input(2381)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Projected CRS: EPSG:2381>\n",
       "Name: Xian 1980 / 3-degree Gauss-Kruger CM 108E\n",
       "Axis Info [cartesian]:\n",
       "- X[north]: Northing (metre)\n",
       "- Y[east]: Easting (metre)\n",
       "Area of Use:\n",
       "- name: China - 106.5°E to 109.5°E onshore\n",
       "- bounds: (106.5, 18.19, 109.5, 42.47)\n",
       "Coordinate Operation:\n",
       "- name: 3-degree Gauss-Kruger CM 108E\n",
       "- method: Transverse Mercator\n",
       "Datum: Xian 1980\n",
       "- Ellipsoid: IAG 1975\n",
       "- Prime Meridian: Greenwich"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyproj.CRS.from_user_input('EPSG:2381')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'+proj=tmerc +lat_0=0 +lon_0=108 +k=1 +x_0=500000 +y_0=0 +ellps=IAU76 +units=m +no_defs +type=crs'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pyproj.CRS.from_user_input(2381).to_proj4()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 设置CRS"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "139    MULTIPOLYGON (((109.47521 18.19770, 108.65521 ...\n",
       "140    POLYGON ((121.77782 24.39427, 121.17563 22.790...\n",
       "Name: geometry, dtype: geometry"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import geopandas as gpd\n",
    "\n",
    "world = gpd.read_file(gpd.datasets.get_path('naturalearth_lowres'))\n",
    "# 利用name字段选择中国区域\n",
    "china = world.loc[world['name'].isin(['China', 'Taiwan']), 'geometry']\n",
    "china"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'init': 'epsg:4326'}"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "china.crs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8defbbac26af4b228b13948343023207",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x25971ac0c08>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%matplotlib widget\n",
    "\n",
    "china.plot(color='red', alpha=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a675501aa55e49d1be0d08007692f59a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.axes._subplots.AxesSubplot at 0x25972e69fc8>"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "china.to_crs(crs='EPSG:2381').plot(color='red', alpha=0.8)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "158ce6b7af4b42ae9c63298b4786381a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(array([-3000000., -2000000., -1000000.,        0.,  1000000.,  2000000.,\n",
       "         3000000.]), <a list of 7 Text xticklabel objects>)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from shapely import geometry\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "cq = gpd.GeoSeries([geometry.Point([106.561203, 29.558078])],\n",
    "              crs='EPSG:4326')\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "china.to_crs(crs='EPSG:2381').plot(ax=ax, color='red', alpha=0.8)\n",
    "cq.plot(ax=ax, color='orange', markersize=100, marker='x')\n",
    "plt.xticks(rotation=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a34ae00ecca943049d0ef4b1d679e84b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "(array([-3000000., -2000000., -1000000.,        0.,  1000000.,  2000000.,\n",
       "         3000000.]), <a list of 7 Text xticklabel objects>)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "fig, ax = plt.subplots()\n",
    "china.to_crs(crs='EPSG:2381').plot(ax=ax, color='red', alpha=0.8)\n",
    "# 先再投影到EPSG:2381\n",
    "cq.to_crs(crs='EPSG:2381').plot(ax=ax, color='orange', markersize=100, marker='x')\n",
    "plt.xticks(rotation=20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9802084047062.387"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 面积单位：平方米\n",
    "china.to_crs(crs='EPSG:2380').area.sum()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "957.6785537419752"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# 长度单位：度，因此计算出的面积无有意义的单位\n",
    "china.area.sum()"
   ]
  }
 ],
 "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}