{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import folium\n",
"import matplotlib\n",
"import numpy as np\n",
"import branca\n",
"from matplotlib.pyplot import imread\n",
"from matplotlib.colors import Normalize\n",
"from matplotlib.colors import ListedColormap\n",
"from folium import raster_layers"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"image = '/notebooks/resources/gpm/gpm_1d.20190530.tif'\n",
"shapefile = '/notebooks/resources/centro-oeste.geojson'"
]
},
{
"cell_type": "code",
"execution_count": 3,
"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"
]
}
],
"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",
"\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)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# create custom colormap\n",
"\n",
"colormap = branca.colormap.StepColormap(\n",
" ['#2b83ba', '#64abb0', '#9dd3a7', '#c7e9ad', '#edf8b9', '#ffedaa', '#fec980', '#f99e59', '#e85b3a', '#d7191c'],\n",
" vmin=0, vmax=10000,\n",
" index=[0, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000],\n",
" caption='rainfall (mm)'\n",
")\n",
"\n",
"colormap"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# recreate custom colormap setting transparency on the first color\n",
"\n",
"change = list()\n",
"alphacolor = list()\n",
"for item in colormap.colors:\n",
" change.append(list(item))\n",
"\n",
"change[0][3] = 0\n",
"\n",
"for item in change:\n",
" alphacolor.append(tuple(item))"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"# set colormap\n",
"colormap = branca.colormap.StepColormap(\n",
" alphacolor\n",
" , vmin=0, vmax=10000\n",
" , index=[0, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000]\n",
" , caption='rainfall (mm)'\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"# set scale\n",
"colormap_scale = colormap.to_linear().scale(0, 1000).to_step(10)\n",
"colormap_scale.caption = 'rainfall(mm)'"
]
},
{
"cell_type": "code",
"execution_count": 59,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 59,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# draw map\n",
"m = folium.Map(\n",
" location = [-16, -63]\n",
" , zoom_start = 5\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 = 1\n",
" , bounds = [ext[2], ext[0]]\n",
" , mercator_project = True\n",
" , colormap = colormap.rgba_floats_tuple\n",
")\n",
" )\n",
"\n",
"m.add_child(colormap_scale)\n",
"\n",
"from rasterstats import zonal_stats\n",
"\n",
"calcs = zonal_stats(\n",
" shapefile, \n",
" image, \n",
" stats=\"count min mean max median\",\n",
" nodata=-999\n",
")\n",
"\n",
"media = str(\"%.2f\" % calcs[0]['mean'])\n",
"\n",
"text = 'Average Rainfall: ' + media + ' mm'\n",
"\n",
"folium.Marker(\n",
" [-15, -53]\n",
" , popup = text\n",
" , tooltip = text\n",
").add_to(m)\n",
"\n",
"from folium import GeoJson\n",
"\n",
"folium.GeoJson(shapefile, name='geojson').add_to(m)\n",
"\n",
"m"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"# m.save('results/map.html')"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"(700, 1500)"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# data.shape"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Driver: GTiff/GeoTIFF\r\n",
"Files: /notebooks/resources/gpm/gpm_1d.20190530.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:05:31 13:35:19\r\n",
" TIFFTAG_DOCUMENTNAME=/NRTPUB/imerg/gis/2019/05/3B-HHR-L.MS.MRG.3IMERG.20190530-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.20190530.tif'"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"79"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# data[0][0]"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"29.003405714285716"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# np.mean(data)"
]
},
{
"cell_type": "code",
"execution_count": 56,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"'0.80'"
]
},
"execution_count": 56,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# from rasterstats import zonal_stats\n",
"# calcs = zonal_stats(\n",
"# shapefile, \n",
"# image, \n",
"# stats=\"count min mean max median\",\n",
"# nodata=-999\n",
"# )\n",
"\n",
"# type(calcs)\n",
"# calcs[0]['mean']\n",
"\n",
"# print(round(calcs[0]['mean']))\n",
"# print(\"%.2f\" % calcs[0]['mean'])\n",
"\n",
"# str(\"%.2f\" % calcs[0]['mean'])"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# from folium import IFrame\n",
"\n",
"# text = 'your text here'\n",
"\n",
"# iframe = folium.IFrame(text, width=700, height=450)\n",
"# popup = folium.Popup(iframe, max_width=3000)\n",
"\n",
"# Text = folium.Marker(location=[40.738, -73.98], popup=popup,\n",
"# icon=folium.Icon(icon_color='green'))\n",
"# m.add_child(Text)"
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
""
],
"text/plain": [
""
]
},
"execution_count": 28,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# from folium import GeoJson\n",
"\n",
"# folium.GeoJson(shapefile, name='geojson').add_to(m)\n",
"\n",
"# m"
]
},
{
"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
}