Browse Source

update

master
yutsuo 7 years ago
parent
commit
1743d66f3d
  1. 157
      Geotiff.ipynb
  2. 0
      HDF5.ipynb
  3. 205
      Untitled1.ipynb
  4. 124
      Untitled2.ipynb

157
Geotiff.ipynb

@ -0,0 +1,157 @@
{
"cells": [
{
"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",
"\n",
"fname = '/notebooks/resources/T22KHG_20190425T132241_TCI_smaller.tif'\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": {
"scrolled": false
},
"outputs": [],
"source": [
"# Second: convert the WKT projection information to a cartopy projection.\n",
"import cartopy.crs as ccrs\n",
"\n",
"projcs = inproj.GetAuthorityCode('PROJCS')\n",
"print('\\n\\n## projcs ##:\\n\\n' + str(projcs))\n",
"\n",
"projection = ccrs.epsg(projcs)\n",
"print('\\n\\n## projection ##:\\n\\n' + str(projection))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# Third (and last): the figure.\n",
"import matplotlib.pyplot as plt\n",
"%matplotlib inline\n",
"\n",
"subplot_kw = dict(projection=projection)\n",
"fig, ax = plt.subplots(figsize=(9, 9), subplot_kw=subplot_kw)\n",
"\n",
"extent = (gt[0], gt[0] + ds.RasterXSize * gt[1],\n",
" gt[3] + ds.RasterYSize * gt[5], gt[3])\n",
"\n",
"img = ax.imshow(data[:3, :, :].transpose((1, 2, 0)), extent=extent, origin='upper')"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import inspect\n",
"\n",
"for i in inspect.getmembers(ds):\n",
" # Ignores anything starting with underscore \n",
" # (that is, private and protected attributes)\n",
" if not i[0].startswith('_'):\n",
" # Ignores methods\n",
" if not inspect.ismethod(i[1]):\n",
" print(i)\n",
"print('\\n')\n",
"# dir(ds)"
]
},
{
"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",
"\n",
"# fname = 'T22KHG_20190425T132241_TCI_smaller.tif'\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(inproj)\n",
"\n",
"# # Second: convert the WKT projection information to a cartopy projection.\n",
"# import cartopy.crs as ccrs\n",
"\n",
"\n",
"\n",
"# projcs = inproj.GetAuthorityCode('PROJCS')\n",
"# print(projcs)\n",
"# projection = ccrs.epsg(projcs)\n",
"# print(projection)\n",
"\n",
"# # Third (and last): the figure.\n",
"# import matplotlib.pyplot as plt\n",
"# %matplotlib inline\n",
"\n",
"# subplot_kw = dict(projection=projection)\n",
"# fig, ax = plt.subplots(figsize=(9, 9), subplot_kw=subplot_kw)\n",
"\n",
"# extent = (gt[0], gt[0] + ds.RasterXSize * gt[1],\n",
"# gt[3] + ds.RasterYSize * gt[5], gt[3])\n",
"\n",
"# img = ax.imshow(data[:3, :, :].transpose((1, 2, 0)), extent=extent, origin='upper')"
]
}
],
"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
}

0
Untitled.ipynb → HDF5.ipynb

205
Untitled1.ipynb

File diff suppressed because one or more lines are too long

124
Untitled2.ipynb

@ -1,124 +0,0 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"ename": "AttributeError",
"evalue": "'NoneType' object has no attribute 'RasterXSize'",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mAttributeError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-1-fe2196e9e159>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 34\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 35\u001b[0m \u001b[0;31m#Get the dimentions of column and row\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 36\u001b[0;31m \u001b[0mnc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mRasterXSize\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 37\u001b[0m \u001b[0mnr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mds\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mRasterYSize\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 38\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mAttributeError\u001b[0m: 'NoneType' object has no attribute 'RasterXSize'"
]
}
],
"source": [
"# pip install --user https://github.com/matplotlib/basemap/archive/master.zip\n",
"\n",
"import os\n",
"from mpl_toolkits.basemap import Basemap\n",
"import matplotlib.pyplot as plt\n",
"from osgeo import gdal\n",
"from osgeo import osr\n",
"import numpy as np\n",
"import math\n",
"\n",
"#Plot setup\n",
"fig= plt.figure(figsize=(10,8))\n",
"\n",
"ax = plt.subplot(111,aspect = 'equal')\n",
"plt.subplots_adjust(left=0.1, bottom=0.1, right=0.9, top=0.9, wspace=0, hspace=0)\n",
"\n",
"#Map setup\n",
"map = Basemap(resolution='f', projection='cyl', llcrnrlon=-75, llcrnrlat=-5, urcrnrlon=-46, urcrnrlat=0)\n",
"\n",
"parallels = np.arange(-50,50,1.)\n",
"meridians = np.arange(0,360,1)\n",
"\n",
"map.drawparallels(parallels,labels=[1,0,0,0],color='w', fontsize=10, fontweight='bold')\n",
"meri = map.drawmeridians(meridians,labels=[0,0,0,1],color='w', fontsize=10, fontweight='bold')\n",
"\n",
"#Load colormap and setup elevation contour levels\n",
"cmap=plt.get_cmap(\"jet\")\n",
"clevs = np.linspace(0, 305, 21)\n",
"\n",
"# Load GeoTiff data\n",
"raster = 'geo_dem_v2_52-50.tif'\n",
"\n",
"ds = gdal.Open(raster)\n",
"\n",
"#Get the dimentions of column and row\n",
"nc = ds.RasterXSize\n",
"nr = ds.RasterYSize\n",
"\n",
"#Read elevation data\n",
"data = ds.ReadAsArray()\n",
"\n",
"#Get Longitude and Latitude info\n",
"geotransform = ds.GetGeoTransform()\n",
"xOrigin = geotransform[0]\n",
"yOrigin = geotransform[3]\n",
"pixelWidth = geotransform[1]\n",
"pixelHeight = geotransform[5]\n",
"\n",
"#Generate Longitude and Latitude array\n",
"lons = xOrigin + np.arange(0, nc)*pixelWidth\n",
"lats = yOrigin + np.arange(0, nr)*pixelHeight\n",
"\n",
"#Contour plot\n",
"x, y = map(*np.meshgrid(lons, lats))\n",
"cs=map.contourf(x, y, data, clevs, cmap=cmap)\n",
"\n",
"map.drawparallels(parallels,labels=[1,0,0,0],color='k', fontsize=10, fontweight='bold')\n",
"meri = map.drawmeridians(meridians,labels=[0,0,0,1],color='k', fontsize=10, fontweight='bold')\n",
"\n",
"cb = map.colorbar(cs, 'bottom', size='5%', pad='10%')\n",
"\n",
"cb.set_label('Elevation (m)', fontsize=12, fontweight='bold')\n",
"cb.ax.tick_params(labelsize=10)\n",
"\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"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
}
Loading…
Cancel
Save