{ "cells": [ { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# pip install colorcet" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import folium\n", "import matplotlib\n", "import colorcet\n", "from matplotlib.pyplot import imread\n", "from matplotlib.colors import Normalize\n", "from folium import raster_layers" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "image = '/notebooks/resources/gpm/gpm_1d.20190531.tif'" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10.0 -180.0\n", "-60.0 -180.0\n", "-60.0 -30.0\n", "10.0 -30.0\n", "ext = [[10.0, -180.0], [-60.0, -180.0], [-60.0, -30.0], [10.0, -30.0]]\n", "geo_ext = [[10.0, -180.0], [-60.0, -180.0], [-60.0, -30.0], [10.0, -30.0]]\n" ] } ], "source": [ "from osgeo import gdal,ogr,osr\n", "\n", "def GetExtent(gt,cols,rows):\n", " ''' Return list of corner coordinates from a geotransform\n", "\n", " @type gt: C{tuple/list}\n", " @param gt: geotransform\n", " @type cols: C{int}\n", " @param cols: number of columns in the dataset\n", " @type rows: C{int}\n", " @param rows: number of rows in the dataset\n", " @rtype: C{[float,...,float]}\n", " @return: coordinates of each corner\n", " '''\n", " ext=[]\n", " xarr=[0,cols]\n", " yarr=[0,rows]\n", "\n", " for px in xarr:\n", " for py in yarr:\n", " x=gt[0]+(px*gt[1])+(py*gt[2])\n", " y=gt[3]+(px*gt[4])+(py*gt[5])\n", " ext.append([y,x])\n", " print (y,x)\n", " yarr.reverse()\n", " return ext\n", "\n", "def ReprojectCoords(coords,src_srs,tgt_srs):\n", " ''' Reproject a list of x,y coordinates.\n", "\n", " @type geom: C{tuple/list}\n", " @param geom: List of [[x,y],...[x,y]] coordinates\n", " @type src_srs: C{osr.SpatialReference}\n", " @param src_srs: OSR SpatialReference object\n", " @type tgt_srs: C{osr.SpatialReference}\n", " @param tgt_srs: OSR SpatialReference object\n", " @rtype: C{tuple/list}\n", " @return: List of transformed [[x,y],...[x,y]] coordinates\n", " '''\n", " trans_coords=[]\n", " transform = osr.CoordinateTransformation( src_srs, tgt_srs)\n", " for x,y in coords:\n", " x,y,z = transform.TransformPoint(x,y)\n", " trans_coords.append([x,y])\n", " return trans_coords\n", "\n", "raster=image\n", "ds=gdal.Open(raster)\n", "\n", "gt=ds.GetGeoTransform()\n", "cols = ds.RasterXSize\n", "rows = ds.RasterYSize\n", "\n", "ext=GetExtent(gt,cols,rows)\n", "print(\"ext = \" + str(ext))\n", "\n", "src_srs=osr.SpatialReference()\n", "src_srs.ImportFromWkt(ds.GetProjection())\n", "# tgt_srs=osr.SpatialReference()\n", "# tgt_srs.ImportFromEPSG(3857)\n", "tgt_srs = src_srs.CloneGeogCS()\n", "\n", "geo_ext=ReprojectCoords(ext,src_srs,tgt_srs)\n", "print(\"geo_ext = \" + str(geo_ext))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: GTiff/GeoTIFF\r\n", "Files: /notebooks/resources/gpm/gpm_1d.20190531.tif\r\n", "Size is 1500, 700\r\n", "Coordinate System is:\r\n", "GEOGCS[\"WGS 84\",\r\n", " DATUM[\"WGS_1984\",\r\n", " SPHEROID[\"WGS 84\",6378137,298.2572326660159,\r\n", " AUTHORITY[\"EPSG\",\"7030\"]],\r\n", " AUTHORITY[\"EPSG\",\"6326\"]],\r\n", " PRIMEM[\"Greenwich\",0],\r\n", " UNIT[\"degree\",0.0174532925199433],\r\n", " AUTHORITY[\"EPSG\",\"4326\"]]\r\n", "Origin = (-180.000000000000000,10.000000000000000)\r\n", "Pixel Size = (0.100000000000000,-0.100000000000000)\r\n", "Metadata:\r\n", " AREA_OR_POINT=Area\r\n", " TIFFTAG_DATETIME=2019:06:01 13:35:13\r\n", " TIFFTAG_DOCUMENTNAME=/NRTPUB/imerg/gis/2019/05/3B-HHR-L.MS.MRG.3IMERG.20190531-S233000-E235959.1410.V06B.1day.tif\r\n", " TIFFTAG_IMAGEDESCRIPTION=DOI=10.5067/GPM/IMERG/3B-HH-L/06 DOIauthority=http://dx.doi.org/ DOIshortName=3IMERGHH_LATE Unit=0.1(mm) ScaleFactor=10\r\n", " TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)\r\n", " TIFFTAG_SOFTWARE=IDL 8.7.2, Harris Geospatial Solutions, Inc.\r\n", " TIFFTAG_XRESOLUTION=100\r\n", " TIFFTAG_YRESOLUTION=100\r\n", "Image Structure Metadata:\r\n", " COMPRESSION=LZW\r\n", " INTERLEAVE=BAND\r\n", "Corner Coordinates:\r\n", "Upper Left (-180.0000000, 10.0000000) (180d 0' 0.00\"W, 10d 0' 0.00\"N)\r\n", "Lower Left (-180.0000000, -60.0000000) (180d 0' 0.00\"W, 60d 0' 0.00\"S)\r\n", "Upper Right ( -30.0000000, 10.0000000) ( 30d 0' 0.00\"W, 10d 0' 0.00\"N)\r\n", "Lower Right ( -30.0000000, -60.0000000) ( 30d 0' 0.00\"W, 60d 0' 0.00\"S)\r\n", "Center (-105.0000000, -25.0000000) (105d 0' 0.00\"W, 25d 0' 0.00\"S)\r\n", "Band 1 Block=1500x2 Type=UInt16, ColorInterp=Gray\r\n" ] } ], "source": [ "!gdalinfo '/notebooks/resources/gpm/gpm_1d.20190531.tif'" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/bin/sh: 1: gdal_edit: not found\r\n" ] } ], "source": [ "!gdal_edit -colorinterp_1 alpha /notebooks/resources/gpm/gpm_1d.20190531.tif" ] }, { "cell_type": "code", "execution_count": 10, "metadata": { "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
" ], "text/plain": [ "" ] }, "execution_count": 10, "metadata": {}, "output_type": "execute_result" } ], "source": [ "m = folium.Map(\n", " location = [-22, -114]\n", " , zoom_start = 2\n", " , control_scale = True \n", " , tiles = 'Stamen Terrain'\n", ")\n", "\n", "data = matplotlib.pyplot.imread(image)\n", "\n", "# Image bounds on the map in the form\n", "# [[lat_min, lon_min], [lat_max, lon_max]]\n", "m.add_child(raster_layers.ImageOverlay(\n", " data\n", " , opacity = 0.7\n", " , bounds = [ext[2], ext[0]]\n", " , mercator_project = True\n", "# , colormap = lambda x: (1, 0, 0, x)\n", " , colormap = colorcet.cm.bmw\n", ")\n", " )\n", "\n", "folium.Marker(\n", " ext[2]\n", " , popup = str(ext[2])\n", " , tooltip = str(ext[2])\n", ").add_to(m)\n", "\n", "folium.Marker(\n", " ext[0]\n", " , popup = str(ext[0])\n", " , tooltip = str(ext[0])\n", ").add_to(m)\n", "\n", "m" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "vars(m)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(data)\n", "print(data.shape)\n", "# data = data.transpose()\n", "# print(data)\n", "# print(data.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# First: read the geotiff image with GDAL.\n", "from osgeo import gdal, osr\n", "\n", "gdal.UseExceptions()\n", "\n", "ds = gdal.Open(image)\n", "data = ds.ReadAsArray()\n", "gt = ds.GetGeoTransform()\n", "proj = ds.GetProjection()\n", "\n", "inproj = osr.SpatialReference()\n", "inproj.ImportFromWkt(proj)\n", "\n", "print('\\n\\n## ds ##:\\n\\n' + str(ds))\n", "print('\\n\\n## data ##:\\n\\n' + str(data))\n", "print('\\n\\n## gt ##:\\n\\n' + str(gt))\n", "print('\\n\\n## proj ##:\\n\\n' + str(proj))\n", "print('\\n\\n## inproj ##:\\n\\n' + str(inproj))" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "q = lambda x: (0, 0, 0, 0)\n", "print(q(9))" ] }, { "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.5.2" } }, "nbformat": 4, "nbformat_minor": 2 }