Creates_mask_from_shapefile.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import shapefile\n",
    "from shapely.geometry import shape, Point"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# creates lat lon array\n",
    "respix = 0.25 #0.1\n",
    "lats = np.arange(1.5, -18.75-respix, -respix)\n",
    "lons = np.arange(278.5, 293.0+respix, respix)-360.0\n",
    "lons, lats = np.meshgrid(lons, lats)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# import shapefile\n",
    "ShapeFile = '/home/joao/Downloads/Cuzco_new2.shp'\n",
    "\n",
    "# creates mask\n",
    "r = shapefile.Reader(ShapeFile)\n",
    "shapes = r.shapes()\n",
    "polygon = shape(shapes[0])\n",
    "Mask = [[1 if polygon.contains(Point(lons[j,i],lats[j,i]))==True else 0 for i in range(lons.shape[1])] for j in range(lons.shape[0])]\n",
    "Mask = np.array(Mask, dtype='int16')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPEAAAD8CAYAAACmVULXAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAATPUlEQVR4nO3df5BdZX3H8ffHTSACUqCxDiaxpJ1UZZyKugIt/eEvNKAVO1NbsP5idBhmwGLHGVH/qNPpP3TsOOoUm+5gShyp1CLUlEmNiFrbUTBBEAkxso0W1lAwSIXCCNndT/84Z/Vms3vPubs3e8+59/OaOcOec0+ec+7o9z7Pec7zPF/ZJiLa6xmDvoGIWJ4EcUTLJYgjWi5BHNFyCeKIlksQR7RcgjhihUjaKulhSfcs8rkkfULSpKS7Jb20TrkJ4oiVcy2wucvn5wGbyu0S4O/qFLqsIJa0WdK+8pfjA8spK2LY2f468JMup1wAfNqF24CTJJ1aVe6qpd6QpDHgauBcYArYJWm77XsX+zfH6Fiv4filXjKiJz/jCZ72U1pOGa975fF+5Ccztc694+6n9gA/6zg0YXuih8utAx7o2J8qjz3Y7R8tOYiBM4FJ2/sBJF1P8UuyaBCv4XjO0quXccmI+m73rcsu45GfzPCtnc+rde7Yqff9zPb4Mi630A9O5bjo5QTxQr8aZ80/SdIlFO171nDcMi4XsfIMzDK7UpebAjZ07K8HDlT9o+U8E9f61bA9YXvc9vhqjl3G5SJWnjGHPFNr64PtwNvLXuqzgZ/a7tqUhuXVxEv61Yhom37VxJI+C7wCWCtpCvgwsBrA9hZgB3A+MAk8CVxcp9zlBPEuYJOkjcCPgAuBtyyjvIjGMWamT9N1bV9U8bmBy3otd8lBbHta0uXATmAM2Gp7z1LLi2iq2eq+pYFaTk2M7R0UTYCIoWRgZpiDOGIUDHVNHDHsDBxq+BJWCeKILozTnI5oNcNMs2M4QRzRTTFiq9kSxBFdiZkFByc2R4I4oouiYytBHNFaxXviBHFEq82mJo5or9TEES1nxEzDl6JLEEdUSHM6osWMeNpjg76NrhLEEV0Ugz3SnI5otXRsRbSYLWacmjii1WYbXhNX/sQslD9G0imSbpF0X/nfk4/ubUYMRtGxtarWNih12gnXcmT+mA8At9reBNxa7kcMnbmOrTrboFReeZH8MRcA28q/twFv6vN9RTTGjFVrG5SltgGeM7eote0HJf3KYicmA0S0WUZsUWSAACYATtQpDV8jIeJIs0PaO/2QpFPLWvhU4OF+3lREUxQTIJodxEu9u+3AO8q/3wF8oT+3E9EsRhzyWK1tUCpr4kXyx1wFfE7Su4D7gTcfzZuMGBSb9g/26JI/JomGYwSo8YM9MmIrogszBDVxxKhresdWgjiiC6MsChDRZsWStc0Ok2bfXcTAZfH4iFYzwztiK2JkNL0mbvZPTMSA2WLWz6i11SFps6R9kiYlHTGFV9IvSfpXSd+RtEfSxVVlpiaO6KLo2OrPkEpJY8DVwLnAFLBL0nbb93acdhlwr+0/kPRsYJ+k62w/vVi5CeKIrvq6xtaZwKTt/QCSrqeYm98ZxAaeJUnACRRz+ae7FZogjuii6Niq/Uy8VtLujv2JcirunHXAAx37U8BZ88r4W4oJRgeAZwF/YrtriuQEcUSFHkZsHbQ93uXzhX4N5s+xfx1wF/Aq4NeBWyT9h+3HFis0HVsRXcyN2Kqz1TAFbOjYX09R43a6GLjRhUngB8ALuhWaII6o0MeF8nYBmyRtlHQMcCFF07nT/ZQzBCU9B3g+sL9boWlOR3Rhw6HZ/tR1tqclXQ7sBMaArbb3SLq0/HwL8FfAtZK+S9H8vtL2wW7lJogjuiia0/1rsNreAeyYd2xLx98HgNf2UmaCOKJC00dsJYgjuujxFdNA1EnjskHSVyXtLYeBXVEeTyqXGAH9HXZ5NNS58jTwPtsvBM4GLpN0OknlEiNitlxnq2oblDoL5T0IzGV7eFzSXoqRJxdQrIIJRSqXrwFXHpW7jBiQond6cMvR1tHTM7Gk04CXALdTM5VL0rhEmw3V8jySTgA+D7zX9mPF+OxqSeMSbTcUS9ZKWk0RwNfZvrE8nFQuMfSGpXdawKeAvbY/2vFRUrnESGh673Sdmvgc4G3AdyXdVR77EEnlEiPAFtNtX2PL9n+y8BQqSCqXGAFNb05nxFZEF214Jk4QR1RIEEe02FC9J44YVUPxnjhiVNkw3adFAY6WBHFEhTSnI1osz8QRQ8AJ4oh2S8dWRIvZeSaOaDkxk97piHbLM3FEi2XsdETbuXgubrIEcUSF9E5HtJjTsRXRfk1vTtdZY2uNpG9J+k6ZAeIvy+PJABEjwVatbVDqtBOeAl5l+8XAGcBmSWeTDBAxAuwhCOIyY/n/lbury80UGSC2lce3AW86KncYMWCzVq1tUGo9sUsaK1e6fBi4xfYRGSCARTNASNotafchnurXfUesGLveNii1OrZszwBnSDoJuEnSi+peIBkgos2MmG1473RPd2f7fykSp22mzAABkAwQMcxccxuUOr3Tzy5rYCQ9E3gN8D2SASJGQZ87tiRtlrRP0qSkBTuDJb1C0l3l26B/ryqzTnP6VGCbpDGKoP+c7ZslfZNkgIhR0Kdqtoyhq4FzgSlgl6Tttu/tOOck4JPAZtv3L5ZttFOdDBB3U6QznX/8EZIBIkZAH18fnQlM2t4PIOl6irc893ac8xbgRtv3F9d25WNqs5/YIwbMwOysam3A2rk3MeV2ybzi1gEPdOxPlcc6/QZwsqSvSbpD0tur7jHDLiO6MVC/Jj5oe7zL5wsVNL+xvgp4GUUr95nANyXdZvv7ixWaII6o0Md3wFPAho799cCBBc45aPsJ4AlJXwdeDCwaxGlOR1Tp3zumXcAmSRslHQNcSPGWp9MXgN+VtErSccBZwN5uhaYmjuiqf+OibU9LuhzYCYwBW23vkXRp+fkW23slfRG4G5gFrrF9T7dyE8QRVfo4ksP2DmDHvGNb5u1/BPhI3TITxBHdGDyblT0iWi5BHNFuDZ+2kyCOqJIgjmix3gZ7DESCOKJC0xfKSxBHVEnvdES7KTVxRIsNetmOGhLEEV0pHVsRrZeaOKLlZgd9A93VnopYrj19p6Sby/2kcYnhN/eeuM42IL3MJ76Cw+c1Jo1LjAS53jYodTNArAdeD1zTcThpXGI0NHzh6bo18ceA93P400HSuEQ0QJ3F498APGz7jqVcwPaE7XHb46s5dilFRAxU05vTdXqnzwHeKOl8YA1woqTPUKZxsf1g0rjE0DKNH3ZZJ7XpB22vt30axcJeX7H9VpLGJUZFw5+Jl/Oe+CqSxiVGwFCNnbb9NYqsiEnjEqNjmII4YiQliCPaa9A9z3UkiCOqNLx3OkEcUSE1cUTbJYgjWizPxBFDIEEc0W4alkUBIqKZUhNHVElzOqLF0rEVMQQSxBEtlyCOaC+R3umIdqu5NE/d52ZJmyXtkzQpadEVYiW9XNKMpD+qKjNBHFGlTyt7SBoDrgbOA04HLpJ0+iLn/TWws87tJYgjqvRveZ4zgUnb+20/DVxPsfTzfO8BPk/NdetqPRNL+iHwODADTNsel3QK8E/AacAPgT+2/Wid8iLapIdXTGsl7e7Yn7A90bG/DnigY38KOOuwa0nrgD8EXgW8vM5Fe+nYeqXtgx37cxkgrirb9h8AruyhvIh2qB/EB22Pd/l8oYnJ80v/GHCl7Rmp3jzm5fROXwC8ovx7G8XaWwniGC7ua+/0FLChY389cGDeOePA9WUArwXOlzRt+18WK7RuEBv4kiQDf182EQ7LACFp0QwQwCUAaziu5uUiGqR/74l3AZskbQR+RLEE9FsOu5S9ce5vSdcCN3cLYKgfxOfYPlAG6i2Svlf3rsuAnwA4Uac0/LV5xJH6NezS9rSkyyl6nceArbb3SLq0/HzLUsqtFcS2D5T/fVjSTRS9bMkA0VA7D9z1879f99wzBngnQ6KPVY/tHcCOeccWDF7b76xTZp1cTMdLetbc38BrgXtIBogYBXVfLzU8A8RzgJvKB+1VwD/a/qKkXSQDRKN01sDdjqV2rk8MwSwm2/uBFy9wPBkgYiS0PogjRl6COI6mhZrL0WcJ4ogWy8oeEUMgQRxHU2dPc5rWR0fTFwVIEEdUSHM6Gicjunow4IEcdSSII6okiCPaayhGbEXzpANrZWm22VGcII7oJs/EEe2X5nT0XT/fDaenuoYEcUS7pSaO1qiq1Ue2pk4QR7RYf1e7PCoSxBFd5D1xHHUrOQFiZJvbbnYU18rFJOkkSTdI+p6kvZJ+S9Ipkm6RdF/535OP9s1GDEI/syIeDXUTqn0c+KLtF1Cst7WXX6Rx2QTcWu5HDJdhWO1S0onA7wHvBCizuT0tKWlcGqaqOZvhmkvT9I6tOjXxrwE/Bv5B0p2SrinXnz4sjQuwaBoXSbsl7T7EU3278YiVotl626DU6dhaBbwUeI/t2yV9nB6azknj0hwL1dTLrZ2HtjNrjhmKjq0pYMr27eX+DRRB/VCZvoWkcYlh1vqOLdv/Azwg6fnloVcD95I0LjEq2t6xVXoPcJ2kY4D9wMUUPwBJ4zKChr4J3WFoBnvYvosi+fF8SeMSw83OogARrdfsGE4QR1QZiuZ0xMgykOZ0NE1GbvWo2TFce+x0xMjq53tiSZsl7ZM0KemIQVOS/lTS3eX2DUlH5AafLzVxRIV+9U5LGgOuBs6lGES1S9J22/d2nPYD4PdtPyrpPIrRjmd1KzdBPOKSkK1CfwdynAlM2t4PIOl64AKKwVPF5exvdJx/G7C+qtAEcUQXxWCP2lG8VtLujv2Jcu7AnHXAAx37U3SvZd8F/FvVRRPEI6hqxNVitfMojdQ6TP0ZSgdtLzQoao4WOLbgL4SkV1IE8e9UXTRBHFGhh5q4yhSwoWN/PXDgiOtJvwlcA5xn+5GqQtM7HdFNf1f22AVskrSxnIdwIcVEop+T9DzgRuBttr9fp9DUxPFzCzWXF2paj1azun9jp21PS7oc2AmMAVtt75F0afn5FuAvgF8GPikJYLqiiZ4gjqjUx0UBbO8Adsw7tqXj73cD7+6lzARxRDdZPD6GyWg1ozs0fHmeBHFElWbHcII4oopmm92ernzFJOn5ku7q2B6T9N5kgIiRYIrBHnW2AamzUN4+22fYPgN4GfAkcBPJABEjQBi53jYovQ72eDXwX7b/m2Lg9rby+DbgTf28sYjGsOttA9LrM/GFwGfLvw/LACFp0QwQwCUAazhuqfcZMTgN752uXROXw8TeCPxzLxewPWF73Pb4ao7t9f4iBqsFz8S91MTnAd+2/VC5/5CkU8taOBkgYmi1vne6w0X8oikNyQARI6Hm83DTO7YkHUexpMiNHYevAs6VdF/52VX9v72IAZtLqNbgIK6bAeJJipkVncceIRkgYhQ0uzWdEVsRVQb5DriOBHFElQRxRIvZMNPs9nSCOKJKauKIlksQR7RYEqpFtJ3BeSaOaC+Tjq2I1sszcUTLJYgj2myw46LrSBBHdGOg4VMRE8QRVVITR7RZhl1GtJvBeU8c0XIZsRXRcnkmjmgxu/G903XX2PpzSXsk3SPps5LWJI1LjIyGr7FVJxfTOuDPgHHbL6LIcH4hSeMSI8F4ZqbWNih1l6xdBTxT0irgOOAASeMSo2BuKmKdbUDqJFT7EfA3wP3Ag8BPbX+JeWlcgEXTuEjaLWn3IZ7q351HrBTP1tsGpE5z+mSKWncj8FzgeElvrXuBpHGJNjPgWdfa6pC0WdI+SZOSjngEVeET5ed3S3ppVZl1mtOvAX5g+8e2D1EsIP/blGlcygsnjUsMJ7tvNbGkMeBqipRIpwMXSTp93mnnAZvK7RLg76rKrRPE9wNnSzpOkigWjN9L0rjEiOhjx9aZwKTt/bafBq6naOV2ugD4tAu3ASfNVZaLqXxPbPt2STcA3wamgTuBCeAE4HOS3kUR6G+uKutxHj34Zd/wBHCw6tyWWUu+UxP96nILeJxHd37ZN6ytefoaSbs79idsT3TsrwMe6NifAs6aV8ZC56yj6I9aUN00Lh8GPjzv8FP0mMbF9rMl7bY93su/a7p8p+Fle3Mfi9NCl1jCOYfpJStiRCzPFLChY389xevaXs85TII4YuXsAjZJ2ijpGIpBU9vnnbMdeHvZS302xSvdRZvSMJix0xPVp7ROvlNUsj0t6XJgJ8XIx62290i6tPx8C7ADOB+YBJ4ELq4qV274DI2I6C7N6YiWSxBHtNyKBnHVkLM2kLRB0lcl7S2nZ15RHm/11ExJY5LulHRzud/q7zNKViyIaw45a4Np4H22XwicDVxWfo+2T828gmIk3py2f5+RsZI1cZ0hZ41n+0Hb3y7/fpzi//jraPHUTEnrgdcD13Qcbu33GTUrGcSLDSdrLUmnAS8Bbqfm1MyG+hjwfqBzFH+bv89IWckg7nk4WZNJOgH4PPBe248N+n6WStIbgIdt3zHoe4mlWcnBHj0PJ2sqSaspAvg62zeWhx+SdKrtB1s2NfMc4I2SzgfWACdK+gzt/T4jZyVr4jpDzhqvnI75KWCv7Y92fNTKqZm2P2h7ve3TKP43+Yrtt9LS7zOKVqwmXmzI2Updv4/OAd4GfFfSXeWxDwFX0ePUzIYbtu8ztDLsMqLlMmIrouUSxBEtlyCOaLkEcUTLJYgjWi5BHNFyCeKIlvt/0A0Gz5Rs4D0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "img = plt.imshow(Mask)\n",
    "plt.colorbar(img)\n",
    "plt.show()"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}