Whitelake in general for Khitrov.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# импортируем программные библиотеки\n",
    "from osgeo import gdal\n",
    "import matplotlib.pyplot as plt\n",
    "import matplotlib as mpl\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "#подгрузка вероятностей, землепользования на 1770 и пустого файла для лэнд юза\n",
    "dataset = gdal.Open('Probability/lu1860.asc')\n",
    "dataset0 = gdal.Open('Probability/P0.asc')\n",
    "dataset1 = gdal.Open('Probability/P1.asc')\n",
    "dataset2 = gdal.Open('Probability/P2.asc')\n",
    "dataset3 = gdal.Open('Probability/P3.asc')\n",
    "dataset4 = gdal.Open('Probability/Empty.asc')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "#determining raster size\n",
    "X=dataset0.RasterXSize\n",
    "Y=dataset0.RasterYSize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-99999., -99999., -99999., ..., -99999., -99999., -99999.],\n",
       "       [-99999., -99999., -99999., ..., -99999., -99999., -99999.],\n",
       "       [-99999., -99999., -99999., ..., -99999., -99999., -99999.],\n",
       "       ...,\n",
       "       [-99999., -99999., -99999., ..., -99999., -99999., -99999.],\n",
       "       [-99999., -99999., -99999., ..., -99999., -99999., -99999.],\n",
       "       [-99999., -99999., -99999., ..., -99999., -99999., -99999.]],\n",
       "      dtype=float32)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Матрица чисел\n",
    "dataset.ReadAsArray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "LU1860=dataset.ReadAsArray()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fields = 3958\n",
      "diff1 = 0.04590770837813757\n",
      "one more\n",
      "Fields = 4902\n",
      "diff1 = 0.03865710160066361\n",
      "one more\n",
      "Fields = 5714\n",
      "diff1 = 0.03242035085563304\n",
      "one more\n",
      "Fields = 6378\n",
      "diff1 = 0.027320347783342032\n",
      "one more\n",
      "Fields = 7056\n",
      "diff1 = 0.022112814525791877\n",
      "one more\n",
      "Fields = 7601\n",
      "diff1 = 0.017926818028203624\n",
      "one more\n",
      "Fields = 8038\n",
      "diff1 = 0.014570340102614515\n",
      "one more\n",
      "Fields = 8388\n",
      "diff1 = 0.011882085471135825\n",
      "one more\n",
      "Fields = 8666\n",
      "diff1 = 0.009746843220989884\n",
      "success\n"
     ]
    }
   ],
   "source": [
    "#открытие растров\n",
    "P0=dataset0.ReadAsArray()\n",
    "P1=dataset1.ReadAsArray()\n",
    "P2=dataset2.ReadAsArray()\n",
    "P3=dataset3.ReadAsArray()\n",
    "LU=dataset4.ReadAsArray()\n",
    "#ввод E0, E1, E2, E3\n",
    "iter0=0\n",
    "iter1=0\n",
    "iter2=0\n",
    "iter3=0\n",
    "count=0\n",
    "#ввод эластичностей\n",
    "E0=0.2177815772600983\n",
    "E1=0.43029512473503995\n",
    "E2=0.24595155709342562\n",
    "E3=3.6788990825688074\n",
    "#добавление эластичностей в пикселях с опр. LU (поощряет сохранение в старых границах)\n",
    "for i in range(0, Y):\n",
    "    for j in range(0, X):\n",
    "        if P0[i][j]!=-99999:\n",
    "            if LU[i][j]==0:\n",
    "                P0[i][j]+=E0\n",
    "            if LU[i][j]==1:\n",
    "                P1[i][j]+=E1\n",
    "            if LU[i][j]==2:\n",
    "                P2[i][j]+=E2\n",
    "            if LU[i][j]==3:\n",
    "                P3[i][j]+=E3\n",
    "            else:\n",
    "                continue\n",
    "#несколько итераций\n",
    "#ввод maxcount\n",
    "maxcount=20000\n",
    "while count<maxcount:\n",
    "    #добавление итерационного параметра\n",
    "    for i in range(0, Y):\n",
    "        for j in range(0, X):\n",
    "            if P0[i][j]!=-99999:\n",
    "                P0[i][j]+=iter0\n",
    "                P1[i][j]+=iter1\n",
    "                P2[i][j]+=iter2\n",
    "                P3[i][j]+=iter3\n",
    "            else:\n",
    "                continue\n",
    "        #первичное распределение лэнд юза\n",
    "    for i in range(0, Y):\n",
    "            for j in range(0, X):\n",
    "                if P0[i][j]!=-99999:\n",
    "                    PG=max(P0[i][j],P1[i][j],P2[i][j],P3[i][j])\n",
    "                    if P1[i][j]>=PG:\n",
    "                        LU[i][j]=1\n",
    "                    else:\n",
    "                        LU[i][j]=4\n",
    "                else:\n",
    "                    LU[i][j]=-99999\n",
    "                if LU1860[i][j]==5:\n",
    "                        LU[i][j]=5\n",
    "                else:\n",
    "                    continue\n",
    "    #подсчёт доли пикселей от общего LU\n",
    "    uses=0\n",
    "    fields=0\n",
    "    for i in range(0, Y):\n",
    "        for j in range(0, X):\n",
    "            if P0[i][j]!=-99999:\n",
    "                uses+=1\n",
    "                if LU[i][j]==1:\n",
    "                    fields+=1\n",
    "    print ('Fields =', fields)\n",
    "    demand1=0.07630802789640234\n",
    "    diff1=demand1-fields/uses\n",
    "    #выведем разницы\n",
    "    print ('diff1 =', diff1)\n",
    "    check1=0\n",
    "    if diff1>0.01 or diff1<-0.01: #число из головы (1%)\n",
    "        iter1=diff1\n",
    "    else:\n",
    "        check1=1\n",
    "    #находим подходящее распределение\n",
    "    if check1==1:\n",
    "        print ('success')\n",
    "        break\n",
    "    else:\n",
    "        print ('one more')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "#импорт созданного грида (без привязки!)\n",
    "#координаты добавлены вручную\n",
    "np.savetxt('C:/Users/olegz/Desktop/CLUE/White Lake project/Predicted/Predicted16.asc', LU, fmt='%4.1f', comments='', header='NCOLS %s \\nNROWS %s \\nXLLCORNER 569875.0000000000 \\nYLLCORNER 6120445.0000000000 \\nCELLSIZE 30.000000 \\nNODATA_VALUE -99999.0'%(X, Y))"
   ]
  },
  {
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}