tif2png.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "c0dda200-58b1-4481-a80f-a87c7dbd35d0",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Saved  heightdata/6701_50m_33.png\n",
      "Saved  heightdata/7002_50m_33.png\n",
      "Saved  heightdata/7806_50m_33.png\n",
      "Saved  heightdata/6502_50m_33.png\n",
      "Saved  heightdata/7404_50m_33.png\n",
      "Saved  heightdata/7304_50m_33.png\n",
      "Saved  heightdata/66M1_50m_33.png\n",
      "Saved  heightdata/6900_50m_33.png\n",
      "Saved  heightdata/7607_50m_33.png\n",
      "Saved  heightdata/7906_50m_33.png\n",
      "Saved  heightdata/7102_50m_33.png\n",
      "Saved  heightdata/6601_50m_33.png\n",
      "Saved  heightdata/7707_50m_33.png\n",
      "Saved  heightdata/67M1_50m_33.png\n",
      "Saved  heightdata/7204_50m_33.png\n",
      "Saved  heightdata/7504_50m_33.png\n",
      "Saved  heightdata/6800_50m_33.png\n",
      "Saved  heightdata/7910_50m_33.png\n",
      "Saved  heightdata/6500_50m_33.png\n",
      "Saved  heightdata/7203_50m_33.png\n",
      "Saved  heightdata/7608_50m_33.png\n",
      "Saved  heightdata/7503_50m_33.png\n",
      "Saved  heightdata/7000_50m_33.png\n",
      "Saved  heightdata/6703_50m_33.png\n",
      "Saved  heightdata/7711_50m_33.png\n",
      "Saved  heightdata/7605_50m_33.png\n",
      "Saved  heightdata/7809_50m_33.png\n",
      "Saved  heightdata/6902_50m_33.png\n",
      "Saved  heightdata/6603_50m_33.png\n",
      "Saved  heightdata/7810_50m_33.png\n",
      "Saved  heightdata/7403_50m_33.png\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[30], line 28\u001b[0m\n\u001b[1;32m     25\u001b[0m         \u001b[38;5;28;01mfor\u001b[39;00m x \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(width):\n\u001b[1;32m     26\u001b[0m             value \u001b[38;5;241m=\u001b[39m img\u001b[38;5;241m.\u001b[39mgetpixel( (x, y) )\n\u001b[0;32m---> 28\u001b[0m             new_pixels\u001b[38;5;241m.\u001b[39mappend(\u001b[38;5;28;43mint\u001b[39;49m\u001b[43m(\u001b[49m\u001b[43mvalue\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m/\u001b[39;49m\u001b[43m \u001b[49m\u001b[43mscale_factor\u001b[49m\u001b[43m)\u001b[49m)\n\u001b[1;32m     29\u001b[0m new_img \u001b[38;5;241m=\u001b[39m Image\u001b[38;5;241m.\u001b[39mnew(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mL\u001b[39m\u001b[38;5;124m'\u001b[39m, img\u001b[38;5;241m.\u001b[39msize)\n\u001b[1;32m     30\u001b[0m new_img\u001b[38;5;241m.\u001b[39mputdata(new_pixels)\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "from PIL import Image\n",
    "import os\n",
    "\n",
    "# Load the TIFF file\n",
    "for dirpath, dirnames, filenames in os.walk('/Users/kvakkefly/projects/joachim/open-elevation/data'):\n",
    "    for filename in filenames:\n",
    "        full_path = os.path.join(dirpath, filename)    \n",
    "        if full_path.endswith('tif'):\n",
    "            with Image.open(full_path) as img:\n",
    "                if img.mode == 'F':\n",
    "                    # Convert to 'RGB' for color images or 'L' for grayscale\n",
    "                    img = img.convert('L')  # or img.convert('L') for grayscale\n",
    "\n",
    "                    # Get image dimensions\n",
    "                    width, height = img.size\n",
    "\n",
    "                    # List to hold new pixel values\n",
    "                    new_pixels = []\n",
    "\n",
    "                    # Scale factor\n",
    "                    scale_factor = 60  # Example scale factor\n",
    "\n",
    "                    # Iterate over each pixel\n",
    "                    for y in range(height):\n",
    "                        for x in range(width):\n",
    "                            value = img.getpixel( (x, y) )\n",
    "                            \n",
    "                            new_pixels.append(int(value / scale_factor))\n",
    "                new_img = Image.new('L', img.size)\n",
    "                new_img.putdata(new_pixels)\n",
    "                # new_img.save('path/to/save/rescaled_rgb_image.png')\n",
    "                new_file_name = f'heightdata/{filename.replace(\".tif\", \".png\")}'\n",
    "                new_img.save(new_file_name)\n",
    "                print(\"Saved \", new_file_name)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "93834c2e-8dea-4f2a-87d8-aafe2b0642bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6701_50m_33.tif\n",
      "Saved  heightdata/6701_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7002_50m_33.tif\n",
      "Saved  heightdata/7002_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7806_50m_33.tif\n",
      "Saved  heightdata/7806_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6502_50m_33.tif\n",
      "Saved  heightdata/6502_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7404_50m_33.tif\n",
      "Saved  heightdata/7404_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7304_50m_33.tif\n",
      "Saved  heightdata/7304_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/66M1_50m_33.tif\n",
      "Saved  heightdata/66M1_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6900_50m_33.tif\n",
      "Saved  heightdata/6900_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7607_50m_33.tif\n",
      "Saved  heightdata/7607_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7906_50m_33.tif\n",
      "Saved  heightdata/7906_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7102_50m_33.tif\n",
      "Saved  heightdata/7102_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6601_50m_33.tif\n",
      "Saved  heightdata/6601_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7707_50m_33.tif\n",
      "Saved  heightdata/7707_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/67M1_50m_33.tif\n",
      "Saved  heightdata/67M1_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7204_50m_33.tif\n",
      "Saved  heightdata/7204_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7504_50m_33.tif\n",
      "Saved  heightdata/7504_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6800_50m_33.tif\n",
      "Saved  heightdata/6800_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7910_50m_33.tif\n",
      "Saved  heightdata/7910_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6500_50m_33.tif\n",
      "Saved  heightdata/6500_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7203_50m_33.tif\n",
      "Saved  heightdata/7203_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7608_50m_33.tif\n",
      "Saved  heightdata/7608_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7503_50m_33.tif\n",
      "Saved  heightdata/7503_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7000_50m_33.tif\n",
      "Saved  heightdata/7000_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6703_50m_33.tif\n",
      "Saved  heightdata/6703_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7711_50m_33.tif\n",
      "Saved  heightdata/7711_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7605_50m_33.tif\n",
      "Saved  heightdata/7605_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7809_50m_33.tif\n",
      "Saved  heightdata/7809_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6902_50m_33.tif\n",
      "Saved  heightdata/6902_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6603_50m_33.tif\n",
      "Saved  heightdata/6603_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7810_50m_33.tif\n",
      "Saved  heightdata/7810_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7403_50m_33.tif\n",
      "Saved  heightdata/7403_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7708_50m_33.tif\n",
      "Saved  heightdata/7708_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7303_50m_33.tif\n",
      "Saved  heightdata/7303_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6400_50m_33.tif\n",
      "Saved  heightdata/6400_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7506_50m_33.tif\n",
      "Saved  heightdata/7506_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6802_50m_33.tif\n",
      "Saved  heightdata/6802_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7705_50m_33.tif\n",
      "Saved  heightdata/7705_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7909_50m_33.tif\n",
      "Saved  heightdata/7909_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7907_50m_33.tif\n",
      "Saved  heightdata/7907_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6600_50m_33.tif\n",
      "Saved  heightdata/6600_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7103_50m_33.tif\n",
      "Saved  heightdata/7103_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/69M1_50m_33.tif\n",
      "Saved  heightdata/69M1_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7706_50m_33.tif\n",
      "Saved  heightdata/7706_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6801_50m_33.tif\n",
      "Saved  heightdata/6801_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7505_50m_33.tif\n",
      "Saved  heightdata/7505_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7003_50m_33.tif\n",
      "Saved  heightdata/7003_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6700_50m_33.tif\n",
      "Saved  heightdata/6700_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/68M1_50m_33.tif\n",
      "Saved  heightdata/68M1_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6503_50m_33.tif\n",
      "Saved  heightdata/6503_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7807_50m_33.tif\n",
      "Saved  heightdata/7807_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6901_50m_33.tif\n",
      "Saved  heightdata/6901_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7405_50m_33.tif\n",
      "Saved  heightdata/7405_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7305_50m_33.tif\n",
      "Saved  heightdata/7305_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7606_50m_33.tif\n",
      "Saved  heightdata/7606_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7101_50m_33.tif\n",
      "Saved  heightdata/7101_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6602_50m_33.tif\n",
      "Saved  heightdata/6602_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7709_50m_33.tif\n",
      "Saved  heightdata/7709_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6401_50m_33.tif\n",
      "Saved  heightdata/6401_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7811_50m_33.tif\n",
      "Saved  heightdata/7811_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6803_50m_33.tif\n",
      "Saved  heightdata/6803_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7507_50m_33.tif\n",
      "Saved  heightdata/7507_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7908_50m_33.tif\n",
      "Saved  heightdata/7908_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7004_50m_33.tif\n",
      "Saved  heightdata/7004_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/65M1_50m_33.tif\n",
      "Saved  heightdata/65M1_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7202_50m_33.tif\n",
      "Saved  heightdata/7202_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6501_50m_33.tif\n",
      "Saved  heightdata/6501_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7609_50m_33.tif\n",
      "Saved  heightdata/7609_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7911_50m_33.tif\n",
      "Saved  heightdata/7911_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6702_50m_33.tif\n",
      "Saved  heightdata/6702_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7001_50m_33.tif\n",
      "Saved  heightdata/7001_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7808_50m_33.tif\n",
      "Saved  heightdata/7808_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/64M1_50m_33.tif\n",
      "Saved  heightdata/64M1_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7104_50m_33.tif\n",
      "Saved  heightdata/7104_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7710_50m_33.tif\n",
      "Saved  heightdata/7710_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/7604_50m_33.tif\n",
      "Saved  heightdata/7604_50m_33.png\n",
      "Reducing file size on  /Users/kvakkefly/projects/joachim/open-elevation/data/6903_50m_33.tif\n",
      "Saved  heightdata/6903_50m_33.png\n"
     ]
    }
   ],
   "source": [
    "from tifffile import imread, imsave\n",
    "import numpy as np\n",
    "\n",
    "def scale_image(image, clip_at):\n",
    "    image = image/clip_at\n",
    "    image = np.clip(image, 0, 1)\n",
    "    image = (image * 255).astype(np.uint8)\n",
    "    return image\n",
    "\n",
    "for dirpath, dirnames, filenames in os.walk('/Users/kvakkefly/projects/joachim/open-elevation/data'):\n",
    "    for filename in filenames:\n",
    "        full_path = os.path.join(dirpath, filename)    \n",
    "        if full_path.endswith('tif'):\n",
    "            image = imread(full_path)\n",
    "            image = scale_image(image, clip_at = 60)\n",
    "            print(\"Reducing file size on \", full_path)\n",
    "\n",
    "            new_img = Image.new('L', (image.shape[1], image.shape[0]))\n",
    "            new_img.putdata(image.flatten().tolist())\n",
    "            new_file_name = f'heightdata/{filename.replace(\".tif\", \".png\")}'\n",
    "            new_img.save(new_file_name)\n",
    "            print(\"Saved \", new_file_name)\n",
    "#            imsave('heightdata/'+filename, image, compression ='zlib')\n",
    "            #imsave('heightdata/'+filename, image)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "3130c6f9-0ab9-4322-bcd8-9ba9feaab144",
   "metadata": {},
   "outputs": [
    {
     "ename": "AttributeError",
     "evalue": "'tuple' object has no attribute 'transpose'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAttributeError\u001b[0m                            Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[39], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m \u001b[43mimage\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mshape\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mtranspose\u001b[49m()\n",
      "\u001b[0;31mAttributeError\u001b[0m: 'tuple' object has no attribute 'transpose'"
     ]
    }
   ],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "30403aa4-74cd-4fca-86bf-365c87d41202",
   "metadata": {},
   "outputs": [],
   "source": [
    "from PIL import Image\n",
    "files = []\n",
    "\n",
    "# for dirpath, dirnames, filenames in os.walk('../open-elevation/data/'):\n",
    "for dirpath, dirnames, filenames in os.walk('heightdata/'):\n",
    "    for filename in filenames:\n",
    "        full_path = os.path.join(dirpath, filename)    \n",
    "        \n",
    "        print(full_path)\n",
    "        if full_path.endswith('tfw'):\n",
    "            image_file_name = full_path.replace(\".tfw\", \".png\")\n",
    "            image = Image.open(image_file_name)\n",
    "            with open(full_path, 'r') as f:\n",
    "                lines = f.readlines()\n",
    "                delta_x = float(lines[0])\n",
    "                delta_y = float(lines[3])\n",
    "                x0 = float(lines[4]) - delta_x * 0.5\n",
    "                y0 = float(lines[5]) - delta_y * 0.5\n",
    "                x1 = x0 + delta_x * image.size[1]\n",
    "                y1 = y0 + delta_y * image.size[0]\n",
    "                files.append({\n",
    "                    'filename': image_file_name,\n",
    "                    'pixel_size': (delta_x, delta_y),\n",
    "                    'coords': ( y1,\n",
    "                                y0,\n",
    "                                x0,\n",
    "                                x1,)\n",
    "                })\n",
    "                print(x0, y0, x1, y1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a86d3def-7fd0-4c60-8b59-4d9219c4d31d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "\n",
    "def latlon_to_utm33(latitude, longitude):\n",
    "    # Constants for WGS84 ellipsoid\n",
    "    a = 6378137.0  # semi-major axis in meters\n",
    "    f = 1 / 298.257223563  # flattening\n",
    "    e_sq = 2 * f - f**2  # square of eccentricity\n",
    "    k0 = 0.9996  # scale factor\n",
    "    lon_origin = 15  # central meridian for UTM Zone 33\n",
    "    \n",
    "    # Convert latitude and longitude from degrees to radians\n",
    "    lat_rad = math.radians(latitude)\n",
    "    lon_rad = math.radians(longitude)\n",
    "    \n",
    "    N = a / math.sqrt(1 - e_sq * math.sin(lat_rad)**2)\n",
    "    T = math.tan(lat_rad)**2\n",
    "    C = e_sq / (1 - e_sq) * math.cos(lat_rad)**2\n",
    "    A = math.cos(lat_rad) * (lon_rad - math.radians(lon_origin))\n",
    "    \n",
    "    M = a * ((1 - e_sq / 4 - 3 * e_sq**2 / 64 - 5 * e_sq**3 / 256) * lat_rad - \n",
    "             (3 * e_sq / 8 + 3 * e_sq**2 / 32 + 45 * e_sq**3 / 1024) * math.sin(2 * lat_rad) + \n",
    "             (15 * e_sq**2 / 256 + 45 * e_sq**3 / 1024) * math.sin(4 * lat_rad) - \n",
    "             (35 * e_sq**3 / 3072) * math.sin(6 * lat_rad))\n",
    "    \n",
    "    x = k0 * N * (A + (1 - T + C) * A**3 / 6 + \n",
    "                  (5 - 18 * T + T**2 + 72 * C - 58 * e_sq) * A**5 / 120)\n",
    "    \n",
    "    y = k0 * (M + N * math.tan(lat_rad) * (A**2 / 2 + (5 - T + 9 * C + 4 * C**2) * A**4 / 24 + \n",
    "                                           (61 - 58 * T + T**2 + 600 * C - 330 * e_sq) * A**6 / 720))\n",
    "    \n",
    "    # Adjust for northern and southern hemispheres\n",
    "    if latitude < 0:\n",
    "        y += 10000000  # 10 million meter offset for southern hemisphere\n",
    "    \n",
    "    x += 500000  # 500,000 meter offset for all UTM zones\n",
    "    \n",
    "    return x, y\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5fa591a1-c7b0-43eb-adf7-9fbc725afceb",
   "metadata": {},
   "outputs": [],
   "source": [
    "def find_file(x, y):\n",
    "    for file in files:\n",
    "        ( y0, y1, x0, x1) = file['coords']\n",
    "        if x0 <= x <= x1 and y0 <= y <= y1:\n",
    "            return file\n",
    "    return None\n",
    "\n",
    "def find_height(latitude, longitude):\n",
    "    \"\"\"\n",
    "    Find the height of the terrain at the given longitude and latitude.\n",
    "    :param longitude:\n",
    "    :param latitude:\n",
    "    :return:\n",
    "    \"\"\"\n",
    "    x, y = latlon_to_utm33(latitude, longitude)\n",
    "    print(\"Found utm33\", x,y)\n",
    "    file = find_file(x, y)\n",
    "    image = Image.open(file['filename'])\n",
    "    ( y0, y1, x0, x1) = file['coords']\n",
    "    (delta_x, delta_y) = file['pixel_size']\n",
    "    x_index = int((x - x0) / delta_x)\n",
    "    y_index = int((y1 - y) / abs(delta_y))\n",
    "    print(x_index, y_index, file['filename'])\n",
    "    return image.getpixel( (y_index, x_index) )\n",
    "    \n",
    "find_height(59.004844, 11.063953)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4270a600-f880-4f6d-b452-efb4089d1245",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}