You can not select more than 25 topics
Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.
244 lines
7.0 KiB
244 lines
7.0 KiB
{ |
|
"cells": [ |
|
{ |
|
"cell_type": "code", |
|
"execution_count": null, |
|
"metadata": { |
|
"scrolled": false |
|
}, |
|
"outputs": [], |
|
"source": [ |
|
"import cartopy.crs as ccrs\n", |
|
"# import cartopy.feature as cf\n", |
|
"import matplotlib.pyplot as plt\n", |
|
"# import numpy as np\n", |
|
"\n", |
|
"# import matplotlib.pyplot as plt\n", |
|
"# from matplotlib import cm\n", |
|
"import numpy as np\n", |
|
"import h5py as h5py\n", |
|
"# import cartopy.crs as ccrs\n", |
|
"import cartopy.feature as cfeature\n", |
|
"from cartopy.io.shapereader import Reader\n", |
|
"from cartopy.feature import ShapelyFeature\n", |
|
"\n", |
|
"%matplotlib notebook\n", |
|
"\n", |
|
"projection = ccrs.PlateCarree()\n", |
|
"\n", |
|
"fig = plt.figure(dpi=150)\n", |
|
"\n", |
|
"ax = plt.axes(projection = projection)\n", |
|
"\n", |
|
"ax.coastlines('10m')\n", |
|
"ax.gridlines()\n", |
|
"# ax.add_feature(cf.BORDERS)\n", |
|
"# ax.add_feature(cf.STATES.with_scale('10m'))\n", |
|
"\n", |
|
"# The extent bounds are specified as an array [[x0, y0], [x1, y1]], \n", |
|
"# where x0 is the left side of the extent, y0 is the top, x1 is the right and y1 is the bottom.\n", |
|
"# extent (x0, x1, y0, y1)\n", |
|
"# extent = [-180,180, -90,90] # world\n", |
|
"# extent = [-90, -30, 20, -60] # south america\n", |
|
"# extent = [-74, -31, 5.5, -33] # brazil\n", |
|
"# extent = [-53.5, -45, -11, -20] # brazil\n", |
|
"# extent = [-48, -47, -16, -17] # brazil\n", |
|
"extent = [-44.996, -44.009, -2.805, -1.809] # brazil\n", |
|
"# extent = [-90, -30, 20, -60]\n", |
|
"# extent = [-100, 30, 0, 80]\n", |
|
"# extent = [-87.35, -79.5, 24.1, 30.8]\n", |
|
"\n", |
|
"\n", |
|
"ax.set_extent(extent, crs=ccrs.PlateCarree())\n", |
|
"\n" |
|
] |
|
}, |
|
{ |
|
"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", |
|
"fname = '/notebooks/resources/T23MNT_20190525T132241_TCI_60m.jp2'\n", |
|
"\n", |
|
"ds = gdal.Open(fname)\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": [ |
|
"projcs = inproj.GetAuthorityCode('PROJCS')\n", |
|
"print('\\n\\n## projcs ##:\\n\\n' + str(projcs))\n", |
|
"\n", |
|
"image_projection = ccrs.epsg(projcs)\n", |
|
"print('\\n\\n## image_projection ##:\\n\\n' + str(image_projection))" |
|
] |
|
}, |
|
{ |
|
"cell_type": "code", |
|
"execution_count": null, |
|
"metadata": {}, |
|
"outputs": [], |
|
"source": [ |
|
"img = plt.imread(fname)" |
|
] |
|
}, |
|
{ |
|
"cell_type": "code", |
|
"execution_count": null, |
|
"metadata": {}, |
|
"outputs": [], |
|
"source": [ |
|
"# import cartopy.crs as ccrs\n", |
|
"# import cartopy.feature as cf\n", |
|
"# import matplotlib.pyplot as plt\n", |
|
"# %matplotlib inline\n", |
|
"\n", |
|
"projection = ccrs.PlateCarree()\n", |
|
"\n", |
|
"fig = plt.figure(dpi=100)\n", |
|
"\n", |
|
"ax = plt.axes(projection = projection)\n", |
|
"\n", |
|
"ax.coastlines('10m')\n", |
|
"ax.gridlines()\n", |
|
"\n", |
|
"# The extent bounds are specified as an array [[x0, y0], [x1, y1]], \n", |
|
"# where x0 is the left side of the extent, y0 is the top, x1 is the right and y1 is the bottom.\n", |
|
"extent = [-44.996, -44.009, -2.805, -1.809] # são luíz do maranhão\n", |
|
"\n", |
|
"ax.set_extent(extent, crs=ccrs.PlateCarree())\n", |
|
"img = plt.imread(fname)\n", |
|
"ax.imshow(img, extent=extent, origin='upper', transform=ccrs.PlateCarree())\n", |
|
"\n", |
|
"# plt.show()\n", |
|
"\n", |
|
"print('\\n\\n## get_extent() ##:\\n\\n' + str(ax.get_extent()))" |
|
] |
|
}, |
|
{ |
|
"cell_type": "code", |
|
"execution_count": null, |
|
"metadata": {}, |
|
"outputs": [], |
|
"source": [] |
|
}, |
|
{ |
|
"cell_type": "code", |
|
"execution_count": null, |
|
"metadata": {}, |
|
"outputs": [], |
|
"source": [ |
|
"# import matplotlib.pyplot as plt\n", |
|
"# from matplotlib import cm\n", |
|
"# import numpy as np\n", |
|
"# import h5py as h5py\n", |
|
"# import cartopy.crs as ccrs\n", |
|
"# import cartopy.feature as cfeature\n", |
|
"# from cartopy.io.shapereader import Reader\n", |
|
"# from cartopy.feature import ShapelyFeature\n", |
|
"# %matplotlib inline\n", |
|
"\n", |
|
"hdf5 = '/notebooks/resources/3B-MO.MS.MRG.3IMERG.20190101-S000000-E235959.01.V06A.HDF5'\n", |
|
"dataset = h5py.File(hdf5,'r')\n", |
|
"\n", |
|
"precip = dataset['Grid/precipitation'][:]\n", |
|
"precip = np.transpose(precip[0])\n", |
|
"\n", |
|
"theLats = dataset['Grid/lat'][:]\n", |
|
"theLons = dataset['Grid/lon'][:]\n", |
|
"\n", |
|
"# [-44.996, -44.009, -2.805, -1.809]\n", |
|
"\n", |
|
"condition = ((theLats < -44) & (theLats > -45))\n", |
|
"extractedLats = np.extract(condition, theLats)\n", |
|
"\n", |
|
"# print(theLats[(theLats < -44) & (theLats > -45)])\n", |
|
"# print(theLats[condition])\n", |
|
"# print(extractedLats)\n", |
|
"\n", |
|
"condition = ((theLons < -1.809) & (theLons > -2.805))\n", |
|
"extractedLons = np.extract(condition, theLons)\n", |
|
"\n", |
|
"# print(theLons[(theLons < -1.809) & (theLons > -2.805)])\n", |
|
"# print(theLons[condition])\n", |
|
"# print(extractedLons)\n", |
|
"\n", |
|
"clevs = np.arange(0,1.26,0.125)\n", |
|
"\n", |
|
"x, y = np.float32(np.meshgrid(theLons, theLats))\n", |
|
"# x, y = np.float32(np.meshgrid(extractedLons, extractedLats))\n", |
|
"\n", |
|
"masked_array = np.ma.masked_where(precip < 0,precip)\n", |
|
"\n", |
|
"cmap = 'nipy_spectral'\n", |
|
"\n", |
|
"plt.title('January 2019 Monthly Average Rain Rate')\n", |
|
"\n", |
|
"font = {'weight' : 'bold', 'size' : 3}\n", |
|
"\n", |
|
"plt.rc('font', **font)\n", |
|
"\n", |
|
"plt.xlim(-2.805, -1.809)\n", |
|
"plt.ylim(-45, -44)\n", |
|
"\n", |
|
"# cs = ax.contourf(np.reshape(x, precip), np.reshape(y, precip), precip, clevs, transform=ccrs.PlateCarree(), cmap=cmap)\n", |
|
"cs = ax.contourf(x, y, precip, clevs, transform=ccrs.PlateCarree(), cmap=cmap)\n", |
|
"\n", |
|
"sm = plt.cm.ScalarMappable(cmap=cmap,norm=plt.Normalize(0,1))\n", |
|
"sm._A = []\n", |
|
"plt.colorbar(sm, ax=ax, label='mm/h', shrink=0.5)\n", |
|
"\n", |
|
"# plt.show()" |
|
] |
|
}, |
|
{ |
|
"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 |
|
}
|
|
|