{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"%%capture\n!pip install gdal\n!pip install rasterio\n!pip install geopandas","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!/usr/bin/env python3\n# -*- coding: utf-8 -*-\n\nimport io\nimport os\nimport re\nimport math\nimport time\nimport argparse\nimport itertools\nimport concurrent.futures\nfrom IPython.display import clear_output\n\nfrom PIL import Image\nfrom PIL import TiffImagePlugin\nImage.MAX_IMAGE_PIXELS = None\n\ntry:\n    import httpx\n    SESSION = httpx.Client()\nexcept ImportError:\n    import requests\n    SESSION = requests.Session()\n\n\nSESSION.headers.update({\n    \"Accept\": \"*/*\",\n    \"Accept-Encoding\": \"gzip, deflate\",\n    \"User-Agent\": \"Mozilla/5.0 (Windows NT 10.0; rv:91.0) Gecko/20100101 Firefox/91.0\",\n})\n\nre_coords_split = re.compile('[ ,;]+')\n\n\nEARTH_EQUATORIAL_RADIUS = 6378137.0\n\nDEFAULT_TMS = 'https://tile.openstreetmap.org/{z}/{x}/{y}.png'\n\n\nWKT_3857 = 'PROJCS[\"WGS 84 / Pseudo-Mercator\",GEOGCS[\"WGS 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUTHORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Mercator_1SP\"],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],EXTENSION[\"PROJ4\",\"+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext +no_defs\"],AUTHORITY[\"EPSG\",\"3857\"]]'\n\n\ndef from4326_to3857(lat, lon):\n    xtile = math.radians(lon) * EARTH_EQUATORIAL_RADIUS\n    ytile = math.log(math.tan(math.radians(45 + lat / 2.0))) * EARTH_EQUATORIAL_RADIUS\n    return (xtile, ytile)\n\n\ndef deg2num(lat, lon, zoom):\n    n = 2 ** zoom\n    xtile = ((lon + 180) / 360 * n)\n    ytile = (1 - math.asinh(math.tan(math.radians(lat))) / math.pi) * n / 2\n    return (xtile, ytile)\n\n\ndef is_empty(im):\n    extrema = im.getextrema()\n    if len(extrema) >= 3:\n        if len(extrema) > 3 and extrema[-1] == (0, 0):\n            return True\n        for ext in extrema[:3]:\n            if ext != (0, 0):\n                return False\n        return True\n    else:\n        return extrema[0] == (0, 0)\n\n\ndef paste_tile(bigim, base_size, tile, corner_xy, bbox):\n    if tile is None:\n        return bigim\n    im = Image.open(io.BytesIO(tile))\n    mode = 'RGB' if im.mode == 'RGB' else 'RGBA'\n    size = im.size\n    if bigim is None:\n        base_size[0] = size[0]\n        base_size[1] = size[1]\n        newim = Image.new(mode, (\n            size[0]*(bbox[2]-bbox[0]), size[1]*(bbox[3]-bbox[1])))\n    else:\n        newim = bigim\n\n    dx = abs(corner_xy[0] - bbox[0])\n    dy = abs(corner_xy[1] - bbox[1])\n    xy0 = (size[0]*dx, size[1]*dy)\n    if mode == 'RGB':\n        newim.paste(im, xy0)\n    else:\n        if im.mode != mode:\n            im = im.convert(mode)\n        if not is_empty(im):\n            newim.paste(im, xy0)\n    im.close()\n    return newim\n\n\ndef get_tile(url):\n    retry = 3\n    while 1:\n        try:\n            r = SESSION.get(url, timeout=60)\n            break\n        except Exception:\n            retry -= 1\n            if not retry:\n                raise\n    if r.status_code == 404:\n        return None\n    elif not r.content:\n        return None\n    r.raise_for_status()\n    return r.content\n\n\ndef print_progress(progress, total, done=False):\n    if done:\n        print('Downloaded image %d/%d' % (progress, total))\n\n\ndef download_extent(\n    source, lat0, lon0, lat1, lon1, zoom,\n    progress_callback=print_progress,\n    callback_interval=0.05\n):\n    x0, y0 = deg2num(lat0, lon0, zoom)\n    x1, y1 = deg2num(lat1, lon1, zoom)\n    if x0 > x1:\n        x0, x1 = x1, x0\n    if y0 > y1:\n        y0, y1 = y1, y0\n    corners = tuple(itertools.product(\n        range(math.floor(x0), math.ceil(x1)),\n        range(math.floor(y0), math.ceil(y1))))\n    totalnum = len(corners)\n    futures = {}\n    done_num = 0\n    progress_callback(done_num, totalnum, False)\n    last_done_num = 0\n    last_callback = time.monotonic()\n    cancelled = False\n    with concurrent.futures.ThreadPoolExecutor(5) as executor:\n        for x, y in corners:\n            future = executor.submit(get_tile, source.format(z=zoom, x=x, y=y))\n            futures[future] = (x, y) \n        bbox = (math.floor(x0), math.floor(y0), math.ceil(x1), math.ceil(y1))\n        bigim = None\n        base_size = [256, 256]\n        while futures:\n            done, not_done = concurrent.futures.wait(\n                futures.keys(), timeout=callback_interval,\n                return_when=concurrent.futures.FIRST_COMPLETED\n            )\n            for fut in done:\n                bigim = paste_tile(bigim, base_size, fut.result(), futures[fut], bbox)\n                del futures[fut]\n                done_num += 1\n            if time.monotonic() > last_callback + callback_interval:\n                try:\n                    progress_callback(done_num, totalnum, (done_num > last_done_num))\n                except TaskCancelled:\n                    for fut in futures.keys():\n                        fut.cancel()\n                    futures.clear()\n                    cancelled = True\n                    break\n                last_callback = time.monotonic()\n                last_done_num = done_num\n    if cancelled:\n        raise TaskCancelled()\n    progress_callback(done_num, totalnum, True)\n\n    xfrac = x0 - bbox[0]\n    yfrac = y0 - bbox[1]\n    x2 = round(base_size[0]*xfrac)\n    y2 = round(base_size[1]*yfrac)\n    imgw = round(base_size[0]*(x1-x0))\n    imgh = round(base_size[1]*(y1-y0))\n    retim = bigim.crop((x2, y2, x2+imgw, y2+imgh))\n    if retim.mode == 'RGBA' and retim.getextrema()[3] == (255, 255):\n        retim = retim.convert('RGB')\n    bigim.close()\n    xp0, yp0 = from4326_to3857(lat0, lon0)\n    xp1, yp1 = from4326_to3857(lat1, lon1)\n    pwidth = abs(xp1 - xp0) / retim.size[0]\n    pheight = abs(yp1 - yp0) / retim.size[1]\n    matrix = (min(xp0, xp1), pwidth, 0, max(yp0, yp1), 0, -pheight)\n    return retim, matrix\n\n\ndef generate_tiffinfo(matrix):\n    ifd = TiffImagePlugin.ImageFileDirectory_v2()\n    # GeoKeyDirectoryTag\n    gkdt = [\n        1, 1,\n        0,  # GeoTIFF 1.0\n        0,  # NumberOfKeys\n    ]\n    # KeyID, TIFFTagLocation, KeyCount, ValueOffset\n    geokeys = [\n        # GTModelTypeGeoKey\n        (1024, 0, 1, 1),  # 2D projected coordinate reference system\n        # GTRasterTypeGeoKey\n        (1025, 0, 1, 1),  # PixelIsArea\n        # GTCitationGeoKey\n        (1026, 34737, 25, 0),\n        # GeodeticCitationGeoKey\n        (2049, 34737, 7, 25),\n        # GeogAngularUnitsGeoKey\n        (2054, 0, 1, 9102),  # degree\n        # ProjectedCRSGeoKey\n        (3072, 0, 1, 3857),\n        # ProjLinearUnitsGeoKey\n        (3076, 0, 1, 9001),  # metre\n    ]\n    gkdt[3] = len(geokeys)\n    ifd.tagtype[34735] = 3  # short\n    ifd[34735] = tuple(itertools.chain(gkdt, *geokeys))\n    # GeoDoubleParamsTag\n    ifd.tagtype[34736] = 12  # double\n    # GeoAsciiParamsTag\n    ifd.tagtype[34737] = 1  # byte\n    ifd[34737] = b'WGS 84 / Pseudo-Mercator|WGS 84|\\x00'\n    a, b, c, d, e, f = matrix\n    # ModelPixelScaleTag\n    ifd.tagtype[33550] = 12  # double\n    # ModelTiepointTag\n    ifd.tagtype[33922] = 12  # double\n    # ModelTransformationTag\n    ifd.tagtype[34264] = 12  # double\n    # This matrix tag should not be used\n    # if the ModelTiepointTag and the ModelPixelScaleTag are already defined\n    if c == 0 and e == 0:\n        ifd[33550] = (b, -f, 0.0)\n        ifd[33922] = (0.0, 0.0, 0.0, a, d, 0.0)\n    else:\n        ifd[34264] = (\n            b, c, 0.0, a,\n            e, f, 0.0, d,\n            0.0, 0.0, 0.0, 0.0,\n            0.0, 0.0, 0.0, 1.0\n        )\n    return ifd\n\n\ndef save_image(img, filename, matrix, **params):\n    wld_ext = {\n        '.gif': '.gfw',\n        '.jpg': '.jgw',\n        '.jpeg': '.jgw',\n        '.jp2': '.j2w',\n        '.png': '.pgw',\n        '.tif': '.tfw',\n        '.tiff': '.tfw',\n    }\n    basename, ext = os.path.splitext(filename)\n    ext = ext.lower()\n    wld_name = basename + wld_ext.get(ext, '.wld')\n    img_params = params.copy()\n    if ext == '.jpg':\n        img_params['quality'] = 92\n        img_params['optimize'] = True\n    elif ext == '.png':\n        img_params['optimize'] = True\n    elif ext.startswith('.tif'):\n        img_params['compression'] = 'tiff_adobe_deflate'\n        img_params['tiffinfo'] = generate_tiffinfo(matrix)\n    img.save(filename, **img_params)\n    if not ext.startswith('.tif'):\n        with open(wld_name, 'w', encoding='utf-8') as f_wld:\n            a, b, c, d, e, f = matrix\n            f_wld.write('\\n'.join(map(str, (b, e, c, f, a, d, ''))))\n    return img\n\n\ndef save_geotiff_gdal(img, filename, matrix):\n    if 'GDAL_DATA' in os.environ:\n        del os.environ['GDAL_DATA']\n    if 'PROJ_LIB' in os.environ:\n        del os.environ['PROJ_LIB']\n\n    import numpy\n    from osgeo import gdal\n    gdal.UseExceptions()\n\n    imgbands = len(img.getbands())\n    driver = gdal.GetDriverByName('GTiff')\n    gtiff = driver.Create(filename, img.size[0], img.size[1],\n        imgbands, gdal.GDT_Byte,\n        options=['COMPRESS=DEFLATE', 'PREDICTOR=2', 'ZLEVEL=9', 'TILED=YES'])\n    gtiff.SetGeoTransform(matrix)\n    gtiff.SetProjection(WKT_3857)\n    for band in range(imgbands):\n        array = numpy.array(img.getdata(band), dtype='u8')\n        array = array.reshape((img.size[1], img.size[0]))\n        band = gtiff.GetRasterBand(band + 1)\n        band.WriteArray(array)\n    gtiff.FlushCache()\n    return img\n\n\ndef save_image_auto(img, filename, matrix, use_gdal=False, **params):\n    ext = os.path.splitext(filename)[1].lower()\n    if ext in ('.tif', '.tiff') and use_gdal:\n        return save_geotiff_gdal(img, filename, matrix)\n    else:\n        return save_image(img, filename, matrix, **params)\n\n\nclass TaskCancelled(RuntimeError):\n    pass\n\n\ndef parse_extent(s):\n    coords_text = re_coords_split.split(s)\n    return (float(coords_text[1]), float(coords_text[0]),\n            float(coords_text[3]), float(coords_text[2]))\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T18:54:16.626598Z","iopub.execute_input":"2023-08-12T18:54:16.62694Z","iopub.status.idle":"2023-08-12T18:54:16.776397Z","shell.execute_reply.started":"2023-08-12T18:54:16.626904Z","shell.execute_reply":"2023-08-12T18:54:16.775008Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import geopandas as gpd\n\n# Replace 'path/to/your/shapefile.shp' with the actual path to your shapefile\nshapefile_path = \"/kaggle/input/200m-filtred-grid/200m_filtred_200k_w_coords.gpkg\"\ngdf = gpd.read_file(shapefile_path)\n\n# Extract rows 7 to 10 and rows 88 to 100\nrows_to_extract = list(range(142059, 145000)) + list(range(161790, 165000)) + list(range(182291, 185000))\ngdf1 = gdf.iloc[rows_to_extract]","metadata":{"execution":{"iopub.status.busy":"2023-08-12T18:54:16.779031Z","iopub.execute_input":"2023-08-12T18:54:16.779423Z","iopub.status.idle":"2023-08-12T18:54:24.517746Z","shell.execute_reply.started":"2023-08-12T18:54:16.779379Z","shell.execute_reply":"2023-08-12T18:54:24.516336Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"gdf.info()","metadata":{"execution":{"iopub.status.busy":"2023-08-12T18:54:24.519778Z","iopub.execute_input":"2023-08-12T18:54:24.520171Z","iopub.status.idle":"2023-08-12T18:54:24.569942Z","shell.execute_reply.started":"2023-08-12T18:54:24.520127Z","shell.execute_reply":"2023-08-12T18:54:24.568937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n\ndef download_images_from_gdf_old(gdf, zoom, output_folder):\n    os.mkdir(output_folder)\n    iterations_to_free_memory = 1000\n    iteration_counter = 0\n    for idx, row in gdf.iterrows():\n        right, left, top, bottom = row['right'], row['left'], row['top'], row['bottom']\n        # Define the TMS URL with the current rectangle's coordinates\n        lon0,lat0, lon1, lat1 = left, top, right, bottom\n        tms_url = 'https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}'\n        # Call the download_extent function to download and merge tiles\n        img, matrix = download_extent(tms_url, lat0, lon0, lat1, lon1, zoom)\n        output_file = f'{output_folder}/rectangle_{idx}.tif'\n        save_image_auto(img, output_file, matrix)\n        resz_file = f'{output_folder}/rectangle_{idx}_res.tif'\n        #!gdalwarp -ts 512 512 -r cubic output_file resz_file\n        !gdalwarp -ts 1024 1024 -r cubic -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=100 {output_file} {resz_file}\n        os.remove(output_file)\n        clear_output()\n        iteration_counter += 1\n        if iteration_counter >= iterations_to_free_memory:\n            time.sleep(60) \n            iteration_counter = 0\n        \n        time.sleep(0.5)","metadata":{"execution":{"iopub.status.busy":"2023-08-12T18:54:24.570996Z","iopub.execute_input":"2023-08-12T18:54:24.571347Z","iopub.status.idle":"2023-08-12T18:54:24.583315Z","shell.execute_reply.started":"2023-08-12T18:54:24.571317Z","shell.execute_reply":"2023-08-12T18:54:24.582526Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import requests\nfrom requests.exceptions import HTTPError\n\ndef download_images_from_gdf(gdf, zoom, output_folder):\n    os.mkdir(output_folder)\n    for idx, row in gdf.iterrows():\n        right, left, top, bottom = row['right'], row['left'], row['top'], row['bottom']\n        # Define the TMS URL with the current rectangle's coordinates\n        lon0, lat0, lon1, lat1 = left, top, right, bottom\n        tms_url = 'https://mt1.google.com/vt/lyrs=s&x={x}&y={y}&z={z}'\n        try:\n            # Call the download_extent function to download and merge tiles\n            img, matrix = download_extent(tms_url, lat0, lon0, lat1, lon1, zoom)\n        except HTTPError as e:\n            print(f\"HTTP Error: {e}\")\n            print(f\"Skipping rectangle {idx}\")\n            continue\n        except TaskCancelled:\n            print(\"Task was cancelled. Stopping the download process.\")\n            break\n\n        output_file = f'{output_folder}/rectangle_{idx}.tif'\n        save_image_auto(img, output_file, matrix)\n        resz_file = f'{output_folder}/rectangle_{idx}_res.tif'\n        !gdalwarp -ts 1024 1024 -r cubic -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=100 {output_file} {resz_file}\n        os.remove(output_file)\n        clear_output()\n","metadata":{"execution":{"iopub.status.busy":"2023-08-12T18:54:24.584604Z","iopub.execute_input":"2023-08-12T18:54:24.585085Z","iopub.status.idle":"2023-08-12T18:54:24.601123Z","shell.execute_reply.started":"2023-08-12T18:54:24.585044Z","shell.execute_reply":"2023-08-12T18:54:24.599325Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"download_images_from_gdf(gdf1, 20, '/kaggle/working/imagery_tiles')","metadata":{"execution":{"iopub.status.busy":"2023-08-12T18:59:42.408111Z","iopub.execute_input":"2023-08-12T18:59:42.408455Z","iopub.status.idle":"2023-08-12T18:59:55.431184Z","shell.execute_reply.started":"2023-08-12T18:59:42.408402Z","shell.execute_reply":"2023-08-12T18:59:55.43001Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n#!gdal_translate -outsize 512 512 '/kaggle/working/cell_100.tif' 'output_resized_imagedd.tif'","metadata":{"execution":{"iopub.status.busy":"2023-07-21T21:28:58.439523Z","iopub.execute_input":"2023-07-21T21:28:58.439978Z","iopub.status.idle":"2023-07-21T21:28:59.681596Z","shell.execute_reply.started":"2023-07-21T21:28:58.439932Z","shell.execute_reply":"2023-07-21T21:28:59.680195Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!gdal_translate -outsize 512 512 -co \"COMPRESS=JPEG\" -co \"JPEG_QUALITY=99\" /kaggle/working/cell_100.tif output_resized_compressed_image.tif\n","metadata":{"execution":{"iopub.status.busy":"2023-07-21T21:07:10.449335Z","iopub.execute_input":"2023-07-21T21:07:10.449796Z","iopub.status.idle":"2023-07-21T21:07:11.732686Z","shell.execute_reply.started":"2023-07-21T21:07:10.449753Z","shell.execute_reply":"2023-07-21T21:07:11.731604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#!gdalwarp -ts 512 512 -r cubic -co COMPRESS=JPEG -co PHOTOMETRIC=YCBCR -co JPEG_QUALITY=100 /kaggle/working/cell_100.tif /kaggle/working/comp_cell_100_gdal_shh.tif","metadata":{"execution":{"iopub.status.busy":"2023-07-21T20:51:49.754407Z","iopub.execute_input":"2023-07-21T20:51:49.755339Z","iopub.status.idle":"2023-07-21T20:51:51.093495Z","shell.execute_reply.started":"2023-07-21T20:51:49.755277Z","shell.execute_reply":"2023-07-21T20:51:51.092084Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}