{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Use GDAL/OGR for raster/vector operations - osgeo\n", "\n", "---\n", "\n", "[Recorded lecture: 1:28:45 - 2:19:43 ](https://youtu.be/kBRW2Z5jX8M)\n", "\n", "---\n", "\n", "**Setting working directory for the OSGeoLive**\n", "\n", "[GDAL](https://gdal.org/) is a translator library for raster and vector geospatial data formats that is released under an X/MIT style Open Source License by the Open Source Geospatial Foundation. As a library, it presents a single raster abstract data model and single vector abstract data model to the calling application for all supported formats. It also comes with a variety of useful command line utilities for data translation and processing (source https://gdal.org/).\n", "\n", "Open the bash terminal, migrate in the directory, and open the jupter-notebook.ipynb\n", "\n", " cd /media/sf_LVM_shared/my_SE_data/exercise\n", " jupyter-lab gdal_osgeo.ipynb\n", "\n", "We are going to use the jupyter-notebook as an editor for running gdal commands. The beauty of using gdal in the bash environment is that we can combine it with other bash utilities such awk, sed, and other." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Start to use GDAL commands\n", "**Explor the files**\n", "\n", "Change directory and list all the files with the extension .tif" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "/media/sf_LVM_shared/my_SE_data/exercise/geodata/cloud/SA_intra.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/cloud/SA_meanannual.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/dem/GMTED2010.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/dem/SA_all_dep_1km.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/dem/SA_are_1km_msk.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/dem/SA_elevation_mn_GMTED2010_mn_msk.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/dem/SA_elevation_mn_GMTED2010_mn.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/875_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/876_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/877_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/878_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/879_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/880_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/881_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/882_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/883_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/884_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/885_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/886_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/887_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/888_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/889_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/890_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/891_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/892_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/893_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/894_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/895_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/896_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/897_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/898_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/899_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/900_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/901_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/902_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/903_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/904_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/905_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/906_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/907_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/908_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/909_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/910_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/911_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/912_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/913_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/914_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/915_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/916_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/917_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/918_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/919_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/920_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/921_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/922_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/923_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/924_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/925_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/926_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/927_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/928_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/929_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/930_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/931_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/932_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/933_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/934_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/935_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/936_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/937_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/938_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/939_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/940_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/941_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/942_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/glad_ard/943_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/071W_41N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/071W_42N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/072W_40N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/072W_41N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/072W_42N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/073W_40N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/073W_41N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/landsat_ct/073W_42N_med_B1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month10.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month11.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month12.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month1.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month2.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month3.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month4.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month5.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month6.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month7.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month8.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/LST/LST_MOYDmax_month9.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/mask/msk_1km.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/SDM/eBirdSampling_filtered_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/SDM/eBirdSampling_filtered.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/temperature/ug_bio_3.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/vegetation/ETmean08-11_01_msk.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/vegetation/ETmean08-11_crop.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/vegetation/ETmean08-11.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/vegetation/ETstdev08-11.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/vegetation/GPPmean08-11.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/vegetation/GPPstdev08-11.tif\n", "/media/sf_LVM_shared/my_SE_data/exercise/geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif\n" ] } ], "source": [ "!ls /media/sf_LVM_shared/my_SE_data/exercise/geodata/*/*.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Retrive the characteristic of the one tif file " ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: GTiff/GeoTIFF\n", "Files: geodata/cloud/SA_intra.tif\n", "Size is 5880, 8400\n", "Coordinate System is:\n", "GEOGCRS[\"WGS 84\",\n", " DATUM[\"World Geodetic System 1984\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " CS[ellipsoidal,2],\n", " AXIS[\"geodetic latitude (Lat)\",north,\n", " ORDER[1],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " AXIS[\"geodetic longitude (Lon)\",east,\n", " ORDER[2],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " ID[\"EPSG\",4326]]\n", "Data axis to CRS axis mapping: 2,1\n", "Origin = (-83.000000000000000,14.000000000000000)\n", "Pixel Size = (0.008333333333333,-0.008333333333333)\n", "Metadata:\n", " AREA_OR_POINT=Area\n", "Image Structure Metadata:\n", " COMPRESSION=DEFLATE\n", " INTERLEAVE=BAND\n", "Corner Coordinates:\n", "Upper Left ( -83.0000000, 14.0000000) ( 83d 0' 0.00\"W, 14d 0' 0.00\"N)\n", "Lower Left ( -83.0000000, -56.0000000) ( 83d 0' 0.00\"W, 56d 0' 0.00\"S)\n", "Upper Right ( -34.0000000, 14.0000000) ( 34d 0' 0.00\"W, 14d 0' 0.00\"N)\n", "Lower Right ( -34.0000000, -56.0000000) ( 34d 0' 0.00\"W, 56d 0' 0.00\"S)\n", "Center ( -58.5000000, -21.0000000) ( 58d30' 0.00\"W, 21d 0' 0.00\"S)\n", "Band 1 Block=5880x1 Type=UInt16, ColorInterp=Gray\n", " NoData Value=65535\n" ] } ], "source": [ "!gdalinfo geodata/cloud/SA_intra.tif" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: GTiff/GeoTIFF\n", "Files: /home/user/jupyter/notebook_gallery/Rasterio/data/RGB.byte.tif\n", "Size is 791, 718\n", "Coordinate System is:\n", "PROJCRS[\"UTM Zone 18, Northern Hemisphere\",\n", " BASEGEOGCRS[\"Unknown datum based upon the WGS 84 ellipsoid\",\n", " DATUM[\"Not_specified_based_on_WGS_84_spheroid\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1],\n", " ID[\"EPSG\",7030]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433,\n", " ID[\"EPSG\",9122]]]],\n", " CONVERSION[\"Transverse Mercator\",\n", " METHOD[\"Transverse Mercator\",\n", " ID[\"EPSG\",9807]],\n", " PARAMETER[\"Latitude of natural origin\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433],\n", " ID[\"EPSG\",8801]],\n", " PARAMETER[\"Longitude of natural origin\",-75,\n", " ANGLEUNIT[\"degree\",0.0174532925199433],\n", " ID[\"EPSG\",8802]],\n", " PARAMETER[\"Scale factor at natural origin\",0.9996,\n", " SCALEUNIT[\"unity\",1],\n", " ID[\"EPSG\",8805]],\n", " PARAMETER[\"False easting\",500000,\n", " LENGTHUNIT[\"metre\",1],\n", " ID[\"EPSG\",8806]],\n", " PARAMETER[\"False northing\",0,\n", " LENGTHUNIT[\"metre\",1],\n", " ID[\"EPSG\",8807]]],\n", " CS[Cartesian,2],\n", " AXIS[\"easting\",east,\n", " ORDER[1],\n", " LENGTHUNIT[\"metre\",1,\n", " ID[\"EPSG\",9001]]],\n", " AXIS[\"northing\",north,\n", " ORDER[2],\n", " LENGTHUNIT[\"metre\",1,\n", " ID[\"EPSG\",9001]]]]\n", "Data axis to CRS axis mapping: 1,2\n", "Origin = (101985.000000000000000,2826915.000000000000000)\n", "Pixel Size = (300.037926675094809,-300.041782729804993)\n", "Metadata:\n", " AREA_OR_POINT=Area\n", "Image Structure Metadata:\n", " INTERLEAVE=PIXEL\n", "Corner Coordinates:\n", "Upper Left ( 101985.000, 2826915.000) ( 78d57'31.14\"W, 25d30'21.91\"N)\n", "Lower Left ( 101985.000, 2611485.000) ( 78d53'53.28\"W, 23d33'53.97\"N)\n", "Upper Right ( 339315.000, 2826915.000) ( 76d35'57.98\"W, 25d33' 3.15\"N)\n", "Lower Right ( 339315.000, 2611485.000) ( 76d34'29.73\"W, 23d36'21.41\"N)\n", "Center ( 220650.000, 2719200.000) ( 77d45'28.46\"W, 24d33'41.70\"N)\n", "Band 1 Block=791x3 Type=Byte, ColorInterp=Red\n", " Min=0.000 Max=255.000 \n", " Minimum=0.000, Maximum=255.000, Mean=29.948, StdDev=52.341\n", " NoData Value=0\n", " Metadata:\n", " STATISTICS_MAXIMUM=255\n", " STATISTICS_MEAN=29.947726688477\n", " STATISTICS_MINIMUM=0\n", " STATISTICS_STDDEV=52.340921626611\n", "Band 2 Block=791x3 Type=Byte, ColorInterp=Green\n", " Min=0.000 Max=255.000 \n", " Minimum=0.000, Maximum=255.000, Mean=44.516, StdDev=56.934\n", " NoData Value=0\n", " Metadata:\n", " STATISTICS_MAXIMUM=255\n", " STATISTICS_MEAN=44.516147889382\n", " STATISTICS_MINIMUM=0\n", " STATISTICS_STDDEV=56.934318291876\n", "Band 3 Block=791x3 Type=Byte, ColorInterp=Blue\n", " Min=0.000 Max=255.000 \n", " Minimum=0.000, Maximum=255.000, Mean=48.113, StdDev=60.113\n", " NoData Value=0\n", " Metadata:\n", " STATISTICS_MAXIMUM=255\n", " STATISTICS_MEAN=48.113056354743\n", " STATISTICS_MINIMUM=0\n", " STATISTICS_STDDEV=60.112778509941\n" ] } ], "source": [ "! gdalinfo /home/user/jupyter/notebook_gallery/Rasterio/data/RGB.byte.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Visualize the image" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Default software rendering mode (use -h if accelerated video card installed).\n", "Loading tools from /usr/bin/openev/tools/Tool_Export.py\n", "Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py\n", "Loading tools from /usr/bin/openev/tools/Tool_DriverList.py\n" ] } ], "source": [ "! /usr/bin/openev/bin/openev geodata/vegetation/*.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Reply to the following questions:\\\n", "What is the pixel size?\\\n", "How do you get min and max pixel values?" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Understanding data type**\n", "\n", "[Numerical System Decimal vs Binary](http://spatial-ecology.net/docs/source/lectures/20220414_Floating_point.pdf).\n", "[Recorded lecture: 1:28:45 - 1:48:00 ](https://youtu.be/kBRW2Z5jX8M)\n", "\n", "\n", "| | Ranges of GDAL data types | | base^bit | Precision | Image Size |\n", "|-------------------|---------------------------|----------------|----------|-----------|------------|\n", "| GDAL data type | Minimum | Maximum | | | |\n", "| Byte | 0 | 255 | 2^8 | | 39M |\n", "| UInt16 | 0 | 65,535 | 2^16 | | 78M |\n", "| Int16, CInt16 | -32,768 | 32,767 | 2^16 | | 78M |\n", "| UInt32 | 0 | 4,294,967,295 | 2^32 | | 155M |\n", "| Int32, CInt32 | -2,147,483,648 | 2,147,483,647 | 2^32 | | 155M |\n", "| Float32, CFloat32 | -3.4E38 | 3.4E38 | 2^38 | Integer -16777216 16777216 (2^24) | 155M |\n", "| Float64, CFloat64 | -1.79E308 | 1.79E308 | 2^308 | | 309M |\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Understanding NoData Value** [Recorded lecture: 1:47:55 - 2:06:00 ](https://youtu.be/kBRW2Z5jX8M)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: GTiff/GeoTIFF\n", "Files: geodata/vegetation/ETmean08-11.tif\n", "Size is 720, 600\n", "Coordinate System is:\n", "GEOGCRS[\"WGS 84\",\n", " DATUM[\"World Geodetic System 1984\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " CS[ellipsoidal,2],\n", " AXIS[\"geodetic latitude (Lat)\",north,\n", " ORDER[1],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " AXIS[\"geodetic longitude (Lon)\",east,\n", " ORDER[2],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " ID[\"EPSG\",4326]]\n", "Data axis to CRS axis mapping: 2,1\n", "Origin = (29.000000000000000,4.000000004000000)\n", "Pixel Size = (0.008333333340000,-0.008333333340000)\n", "Metadata:\n", " AREA_OR_POINT=Area\n", "Image Structure Metadata:\n", " COMPRESSION=LZW\n", " INTERLEAVE=BAND\n", "Corner Coordinates:\n", "Upper Left ( 29.0000000, 4.0000000) ( 29d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Left ( 29.0000000, -1.0000000) ( 29d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Upper Right ( 35.0000000, 4.0000000) ( 35d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Right ( 35.0000000, -1.0000000) ( 35d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Center ( 32.0000000, 1.5000000) ( 32d 0' 0.00\"E, 1d30' 0.00\"N)\n", "Band 1 Block=720x2 Type=Float32, ColorInterp=Gray\n", " Computed Min/Max=0.547,7.492\n", " NoData Value=-3.39999995214436425e+38\n" ] } ], "source": [ "!gdalinfo -mm geodata/vegetation/ETmean08-11.tif" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "!gdal_edit.py -a_nodata -9999 geodata/vegetation/ETmean08-11.tif" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: GTiff/GeoTIFF\n", "Files: geodata/vegetation/ETmean08-11.tif\n", "Size is 720, 600\n", "Coordinate System is:\n", "GEOGCRS[\"WGS 84\",\n", " DATUM[\"World Geodetic System 1984\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " CS[ellipsoidal,2],\n", " AXIS[\"geodetic latitude (Lat)\",north,\n", " ORDER[1],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " AXIS[\"geodetic longitude (Lon)\",east,\n", " ORDER[2],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " ID[\"EPSG\",4326]]\n", "Data axis to CRS axis mapping: 2,1\n", "Origin = (29.000000000000000,4.000000004000000)\n", "Pixel Size = (0.008333333340000,-0.008333333340000)\n", "Metadata:\n", " AREA_OR_POINT=Area\n", "Image Structure Metadata:\n", " COMPRESSION=LZW\n", " INTERLEAVE=BAND\n", "Corner Coordinates:\n", "Upper Left ( 29.0000000, 4.0000000) ( 29d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Left ( 29.0000000, -1.0000000) ( 29d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Upper Right ( 35.0000000, 4.0000000) ( 35d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Right ( 35.0000000, -1.0000000) ( 35d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Center ( 32.0000000, 1.5000000) ( 32d 0' 0.00\"E, 1d30' 0.00\"N)\n", "Band 1 Block=720x2 Type=Float32, ColorInterp=Gray\n", " Computed Min/Max=-339999995214436424907732413799364296704.000,7.492\n", " NoData Value=-9999\n" ] } ], "source": [ "!gdalinfo -mm geodata/vegetation/ETmean08-11.tif" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [], "source": [ "!gdal_edit.py -a_nodata -339999995214436424907732413799364296704.000 geodata/vegetation/ETmean08-11.tif" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: GTiff/GeoTIFF\n", "Files: geodata/vegetation/ETmean08-11.tif\n", "Size is 720, 600\n", "Coordinate System is:\n", "GEOGCRS[\"WGS 84\",\n", " DATUM[\"World Geodetic System 1984\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " CS[ellipsoidal,2],\n", " AXIS[\"geodetic latitude (Lat)\",north,\n", " ORDER[1],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " AXIS[\"geodetic longitude (Lon)\",east,\n", " ORDER[2],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " ID[\"EPSG\",4326]]\n", "Data axis to CRS axis mapping: 2,1\n", "Origin = (29.000000000000000,4.000000004000000)\n", "Pixel Size = (0.008333333340000,-0.008333333340000)\n", "Metadata:\n", " AREA_OR_POINT=Area\n", "Image Structure Metadata:\n", " COMPRESSION=LZW\n", " INTERLEAVE=BAND\n", "Corner Coordinates:\n", "Upper Left ( 29.0000000, 4.0000000) ( 29d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Left ( 29.0000000, -1.0000000) ( 29d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Upper Right ( 35.0000000, 4.0000000) ( 35d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Right ( 35.0000000, -1.0000000) ( 35d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Center ( 32.0000000, 1.5000000) ( 32d 0' 0.00\"E, 1d30' 0.00\"N)\n", "Band 1 Block=720x2 Type=Float32, ColorInterp=Gray\n", " Computed Min/Max=0.547,7.492\n", " NoData Value=-3.39999995214436425e+38\n" ] } ], "source": [ "!gdalinfo -mm geodata/vegetation/ETmean08-11.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "[Recorded lecture: 2:06:00 - 2:19:43 ](https://youtu.be/kBRW2Z5jX8M) + [Recorded lecture: 42:18 - 1:46:40](https://youtu.be/Yy5UaCOvuYo)\n", "\n", "**Calculate minimum and maximum values for all the images**" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " Computed Min/Max=-10.000,5.000\n", " Computed Min/Max=0.638,6.957\n", " Computed Min/Max=0.547,7.492\n", " Computed Min/Max=0.209,3.332\n", " Computed Min/Max=0.000,9.738\n", " Computed Min/Max=0.000,4.220\n", " Computed Min/Max=0.000,10000.000\n" ] } ], "source": [ "%%bash \n", "for file in geodata/vegetation/*.tif ; do \n", "gdalinfo -mm $file | grep Computed\n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Create a Coefficient of variation**\n", "\n", "GPPmean08-11.tif temporal mean of Gross Primary Productivity\\\n", "GPPstdev08-11.tif temporal standard deviation of Gross Primary Productivity" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.. 0.. 0.. 1.. 1.. 1.. 2.. 2.. 2.. 3.. 3.. 3.. 4.. 4.. 4.. 5.. 5.. 5.. 6.. 6.. 6.. 7.. 7.. 7.. 8.. 8.. 8.. 9.. 9.. 9.. 10.. 10.. 10.. 11.. 11.. 11.. 12.. 12.. 12.. 13.. 13.. 13.. 14.. 14.. 14.. 15.. 15.. 15.. 16.. 16.. 16.. 17.. 17.. 17.. 18.. 18.. 18.. 19.. 19.. 19.. 20.. 20.. 20.. 21.. 21.. 21.. 22.. 22.. 22.. 23.. 23.. 23.. 24.. 24.. 24.. 25.. 25.. 25.. 26.. 26.. 26.. 27.. 27.. 27.. 28.. 28.. 28.. 29.. 29.. 29.. 30.. 30.. 30.. 31.. 31.. 31.. 32.. 32.. 32.. 33.. 33.. 33.. 34.. 34.. 34.. 35.. 35.. 35.. 36.. 36.. 36.. 37.. 37.. 37.. 38.. 38.. 38.. 39.. 39.. 39.. 40.. 40.. 40.. 41.. 41.. 41.. 42.. 42.. 42.. 43.. 43.. 43.. 44.. 44.. 44.. 45.. 45.. 45.. 46.. 46.. 46.. 47.. 47.. 47.. 48.. 48.. 48.. 49.. 49.. 49.. 50.. 50.. 50.. 51.. 51.. 51.. 52.. 52.. 52.. 53.. 53.. 53.. 54.. 54.. 54.. 55.. 55.. 55.. 56.. 56.. 56.. 57.. 57.. 57.. 58.. 58.. 58.. 59.. 59.. 59.. 60.. 60.. 60.. 61.. 61.. 61.. 62.. 62.. 62.. 63.. 63.. 63.. 64.. 64.. 64.. 65.. 65.. 65.. 66.. 66.. 66.. 67.. 67.. 67.. 68.. 68.. 68.. 69.. 69.. 69.. 70.. 70.. 70.. 71.. 71.. 71.. 72.. 72.. 72.. 73.. 73.. 73.. 74.. 74.. 74.. 75.. 75.. 75.. 76.. 76.. 76.. 77.. 77.. 77.. 78.. 78.. 78.. 79.. 79.. 79.. 80.. 80.. 80.. 81.. 81.. 81.. 82.. 82.. 82.. 83.. 83.. 83.. 84.. 84.. 84.. 85.. 85.. 85.. 86.. 86.. 86.. 87.. 87.. 87.. 88.. 88.. 88.. 89.. 89.. 89.. 90.. 90.. 90.. 91.. 91.. 91.. 92.. 92.. 92.. 93.. 93.. 93.. 94.. 94.. 94.. 95.. 95.. 95.. 96.. 96.. 96.. 97.. 97.. 97.. 98.. 98.. 98.. 99.. 99.. 99.. 100 - Done\n" ] } ], "source": [ "%%bash\n", "gdal_calc.py --type=Float32 --overwrite -A geodata/vegetation/GPPstdev08-11.tif -B geodata/vegetation/GPPmean08-11.tif \\\n", "--co=COMPRESS=DEFLATE --co=ZLEVEL=9 --calc=\"( A.astype(float) / ( B.astype(float) + 0.000000001 ) )\" --outfile=geodata/vegetation/GPPcv08-11.tif" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Driver: GTiff/GeoTIFF\n", "Files: geodata/vegetation/GPPcv08-11.tif\n", "Size is 720, 600\n", "Coordinate System is:\n", "GEOGCRS[\"WGS 84\",\n", " DATUM[\"World Geodetic System 1984\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " CS[ellipsoidal,2],\n", " AXIS[\"geodetic latitude (Lat)\",north,\n", " ORDER[1],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " AXIS[\"geodetic longitude (Lon)\",east,\n", " ORDER[2],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " ID[\"EPSG\",4326]]\n", "Data axis to CRS axis mapping: 2,1\n", "Origin = (29.000000000000000,4.000000004000000)\n", "Pixel Size = (0.008333333340000,-0.008333333340000)\n", "Metadata:\n", " AREA_OR_POINT=Area\n", "Image Structure Metadata:\n", " COMPRESSION=DEFLATE\n", " INTERLEAVE=BAND\n", "Corner Coordinates:\n", "Upper Left ( 29.0000000, 4.0000000) ( 29d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Left ( 29.0000000, -1.0000000) ( 29d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Upper Right ( 35.0000000, 4.0000000) ( 35d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Right ( 35.0000000, -1.0000000) ( 35d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Center ( 32.0000000, 1.5000000) ( 32d 0' 0.00\"E, 1d30' 0.00\"N)\n", "Band 1 Block=720x2 Type=Float32, ColorInterp=Gray\n", " Computed Min/Max=0.000,1.476\n", " NoData Value=3.40282346600000016e+38\n" ] } ], "source": [ "! gdalinfo -mm geodata/vegetation/GPPcv08-11.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Change pixel resolution**\n", "\n", "Looping trough the images" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month1.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month1.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month1.tif to destination geodata/LST/LST_MOYDmax_month1_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month2.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month2.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month2.tif to destination geodata/LST/LST_MOYDmax_month2_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month3.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month3.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month3.tif to destination geodata/LST/LST_MOYDmax_month3_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month4.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month4.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month4.tif to destination geodata/LST/LST_MOYDmax_month4_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month5.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month5.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month5.tif to destination geodata/LST/LST_MOYDmax_month5_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month6.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month6.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month6.tif to destination geodata/LST/LST_MOYDmax_month6_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month7.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month7.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month7.tif to destination geodata/LST/LST_MOYDmax_month7_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month8.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month8.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month8.tif to destination geodata/LST/LST_MOYDmax_month8_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month9.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month9.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month9.tif to destination geodata/LST/LST_MOYDmax_month9_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month10.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month10.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month10.tif to destination geodata/LST/LST_MOYDmax_month10_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month11.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month11.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month11.tif to destination geodata/LST/LST_MOYDmax_month11_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n", "Creating output file that is 720P x 600L.\n", "Processing geodata/LST/LST_MOYDmax_month12.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month12.tif.\n", "Copying nodata values from source geodata/LST/LST_MOYDmax_month12.tif to destination geodata/LST/LST_MOYDmax_month12_crop.tif.\n", "...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do \n", " filename=$(basename $file .tif )\n", " gdalwarp -overwrite -te 29 -1 35.0000000048 4.000000004 -tr 0.008333333340000 0.008333333340000 -co COMPRESS=DEFLATE -co ZLEVEL=9 $file geodata/LST/${filename}_crop.tif \n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Get value at one pixel/line image-location**\n", "\n", "Looping trough the images" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 50.8218383789062\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 52.4714660644531\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 50.2210998535156\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 45.1343078613281\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 43.8952941894531\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 45.8667602539062\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 47.9639892578125\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 49.0910339355469\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 49.9898986816406\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 47.5579528808594\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 45.6031494140625\n", "Report:\n", " Location: (1079P,0L)\n", " Band 1:\n", " Value: 47.5201110839844\n" ] } ], "source": [ "%%bash \n", "for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do \n", " gdallocationinfo $file 1079 0 \n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Get value at lat/long image-location**\n", "\n", "Looping trough the images" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "37.4764404296875\n", "40.2584228515625\n", "38.2513122558594\n", "33.0213928222656\n", "30.0732116699219\n", "29.0838928222656\n", "28.9105224609375\n", "29.6008605957031\n", "31.0305480957031\n", "30.5758056640625\n", "30.636962890625\n", "33.5576171875\n" ] } ], "source": [ "%%bash\n", "for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do \n", " gdallocationinfo -valonly -geoloc $file 33.3 2 \n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Get value at multiple lat/long image-location**\\" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "37.4022827148438\n", "35.0345764160156\n", "\n", "40.3694458007812\n", "37.824951171875\n", "\n", "38.5549926757812\n", "36.6663208007812\n", "\n", "32.7738952636719\n", "32.6803283691406\n", "\n", "29.0351257324219\n", "30.0338134765625\n", "\n", "29.025634765625\n", "29.3355102539062\n", "\n", "29.0183715820312\n", "29.4237365722656\n", "\n", "29.1660461425781\n", "29.3502197265625\n", "\n", "29.6674194335938\n", "29.7001647949219\n", "\n", "28.8177490234375\n", "29.066650390625\n", "\n", "29.2975463867188\n", "29.015869140625\n", "\n", "32.8586730957031\n", "31.3289184570312\n", "\n" ] } ], "source": [ "%%bash\n", "# Create the lat long file \n", "echo 32.5 2.5 > geodata/LST/x_y.txt\n", "echo 31.1 2.1 >> geodata/LST/x_y.txt\n", "# looping trough the images\n", "for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do \n", " gdallocationinfo -valonly -geoloc $file < geodata/LST/x_y.txt \n", " echo \"\"\n", "done" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**From Image to text and from txt to image**" ] }, { "cell_type": "code", "execution_count": 24, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Input file size is 720, 600\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "Input file size is 720, 600\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "# from tif to 3 col txt XYZ\n", "gdal_translate -of XYZ geodata/vegetation/ETmean08-11.tif geodata/vegetation/ETmean08-11_crop.txt\n", "# masking\n", "awk '{if ($3>2 && $3<4) {print $1,$2,50 } else {print}}' geodata/vegetation/ETmean08-11_crop.txt > geodata/vegetation/ETmean08-11_crop_trh.txt\n", "# back to tif\n", "gdal_translate -co COMPRESS=DEFLATE -co ZLEVEL=9 -ot Byte geodata/vegetation/ETmean08-11_crop_trh.txt geodata/vegetation/ETmean08-11_crop_trh.tif " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**The use of VRT to stack images and tiling**" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "Driver: VRT/Virtual Raster\n", "Files: geodata/vegetation/stack.vrt\n", " geodata/vegetation/ETmean08-11.tif\n", " geodata/vegetation/ETstdev08-11.tif\n", "Size is 720, 600\n", "Coordinate System is:\n", "GEOGCRS[\"WGS 84\",\n", " DATUM[\"World Geodetic System 1984\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " CS[ellipsoidal,2],\n", " AXIS[\"geodetic latitude (Lat)\",north,\n", " ORDER[1],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " AXIS[\"geodetic longitude (Lon)\",east,\n", " ORDER[2],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " ID[\"EPSG\",4326]]\n", "Data axis to CRS axis mapping: 2,1\n", "Origin = (29.000000000000000,4.000000004000000)\n", "Pixel Size = (0.008333333340000,-0.008333333340000)\n", "Corner Coordinates:\n", "Upper Left ( 29.0000000, 4.0000000) ( 29d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Left ( 29.0000000, -1.0000000) ( 29d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Upper Right ( 35.0000000, 4.0000000) ( 35d 0' 0.00\"E, 4d 0' 0.00\"N)\n", "Lower Right ( 35.0000000, -1.0000000) ( 35d 0' 0.00\"E, 1d 0' 0.00\"S)\n", "Center ( 32.0000000, 1.5000000) ( 32d 0' 0.00\"E, 1d30' 0.00\"N)\n", "Band 1 Block=128x128 Type=Float32, ColorInterp=Undefined\n", " NoData Value=-3.39999995214436387e+38\n", "Band 2 Block=128x128 Type=Float32, ColorInterp=Undefined\n", " NoData Value=-3.39999995214436387e+38\n" ] } ], "source": [ "%%bash\n", "# stck the tif\n", "gdalbuildvrt -overwrite -separate geodata/vegetation/stack.vrt geodata/vegetation/ETmean08-11.tif geodata/vegetation/ETstdev08-11.tif\n", "gdalinfo geodata/vegetation/stack.vrt" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Input file size is 720, 600\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "Input file size is 720, 600\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "Input file size is 720, 600\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "Input file size is 720, 600\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "Input file size is 720, 600\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "# split in tiles \n", "gdal_translate -srcwin 0 0 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_UL.tif\n", "gdal_translate -srcwin 0 300 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_UR.tif\n", "gdal_translate -srcwin 360 0 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_LL.tif\n", "gdal_translate -srcwin 360 300 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_LR.tif\n", " \n", "# recompose the image \n", "gdalbuildvrt -overwrite geodata/vegetation/ETmosaic.vrt geodata/vegetation/stack_UL.tif geodata/vegetation/stack_UR.tif geodata/vegetation/stack_LL.tif geodata/vegetation/stack_LR.tif\n", "gdal_translate -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/vegetation/ETmosaic.vrt geodata/vegetation/ETmosaic.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Create shp border tiles**" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Creating new index file...\n" ] } ], "source": [ "%%bash\n", "rm -f geodata/vegetation/tiles.*\n", "gdaltindex geodata/vegetation/tiles.shp geodata/vegetation/stack_??.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Surface distance/buffer**" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "# Continues distance surface\n", "gdal_proximity.py -co COMPRESS=DEFLATE -co ZLEVEL=9 -values 0 -distunits PIXEL -ot UInt32 geodata/vegetation/ETmean08-11_crop_trh.tif geodata/vegetation/ETmean08-11_crop_proximity.tif" ] }, { "cell_type": "code", "execution_count": 29, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "# Fix buffer surface\n", "gdal_proximity.py -fixed-buf-val 4 -maxdist 4 -nodata 10 -co COMPRESS=DEFLATE -co ZLEVEL=9 -values 0 -distunits PIXEL -ot Byte geodata/vegetation/ETmean08-11_crop_trh.tif geodata/vegetation/ETmean08-11_crop_buffer.tif" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Topography variables**" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n", "0...10...20...30...40...50...60...70...80...90...100 - done.\n" ] } ], "source": [ "%%bash\n", "# calculate slope \n", "gdaldem slope -s 111120 -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/slope.tif \n", " \n", "# calculate apect\n", "gdaldem aspect -zero_for_flat -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/aspect.tif \n", " \n", "# calculate Terrain Ruggedness Index TRI \n", "gdaldem TRI -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/tri.tif \n", " \n", "# calculate Topographic Position Index TPI\n", "gdaldem TPI -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/tpi.tif \n", " \n", "# calculate Roughness \n", "gdaldem roughness -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/roughness.tif " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Start to use OGR Commands\n", "**Select a polygon and create a new shape file**" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INFO: Open of `geodata/shp/TM_WORLD_BORDERS.shp'\n", " using driver `ESRI Shapefile' successful.\n", "\n", "Layer name: TM_WORLD_BORDERS\n", "Metadata:\n", " DBF_DATE_LAST_UPDATE=2021-02-09\n", "Geometry: Polygon\n", "Feature Count: 246\n", "Extent: (-180.000000, -90.000000) - (180.000000, 83.623596)\n", "Layer SRS WKT:\n", "GEOGCRS[\"WGS 84\",\n", " DATUM[\"World Geodetic System 1984\",\n", " ELLIPSOID[\"WGS 84\",6378137,298.257223563,\n", " LENGTHUNIT[\"metre\",1]]],\n", " PRIMEM[\"Greenwich\",0,\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " CS[ellipsoidal,2],\n", " AXIS[\"latitude\",north,\n", " ORDER[1],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " AXIS[\"longitude\",east,\n", " ORDER[2],\n", " ANGLEUNIT[\"degree\",0.0174532925199433]],\n", " ID[\"EPSG\",4326]]\n", "Data axis to CRS axis mapping: 2,1\n", "FIPS: String (2.0)\n", "ISO2: String (2.0)\n", "ISO3: String (3.0)\n", "UN: Integer (3.0)\n", "NAME: String (50.0)\n", "AREA: Integer (7.0)\n", "POP2005: Integer64 (10.0)\n", "REGION: Integer (3.0)\n", "SUBREGION: Integer (3.0)\n", "LON: Real (8.3)\n", "LAT: Real (7.3)\n", "OGRFeature(TM_WORLD_BORDERS):0\n", " FIPS (String) = AC\n", " ISO2 (String) = AG\n", " ISO3 (String) = ATG\n", " UN (Integer) = 28\n", " NAME (String) = Antigua and Barbuda\n", " AREA (Integer) = 44\n", " POP2005 (Integer64) = 83039\n", " REGION (Integer) = 19\n", " SUBREGION (Integer) = 29\n", " LON (Real) = -61.783\n", " LAT (Real) = 17.078\n", "\n", "OGRFeature(TM_WORLD_BORDERS):1\n", " FIPS (String) = AG\n", " ISO2 (String) = DZ\n", " ISO3 (String) = DZA\n", " UN (Integer) = 12\n", " NAME (String) = Algeria\n", " AREA (Integer) = 238174\n", " POP2005 (Integer64) = 32854159\n", " REGION (Integer) = 2\n", " SUBREGION (Integer) = 15\n", " LON (Real) = 2.632\n", " LAT (Real) = 28.163\n", "\n", "OGRFeature(TM_WORLD_BORDERS):2\n", " FIPS (String) = AJ\n", " ISO2 (String) = AZ\n", " ISO3 (String) = AZE\n", " UN (Integer) = 31\n", " NAME (String) = Azerbaijan\n", " AREA (Integer) = 8260\n", " POP2005 (Integer64) = 8352021\n", " REGION (Integer) = 142\n", " SUBREGION (Integer) = 145\n", " LON (Real) = 47.395\n", " LAT (Real) = 40.430\n", "\n", "OGRFeature(TM_WORLD_BORDERS):3\n", " FIPS (String) = AL\n", " ISO2 (String) = AL\n", " ISO3 (String) = ALB\n", " UN (Integer) = 8\n" ] } ], "source": [ "%%bash\n", "ogrinfo -al -geom=NO geodata/shp/TM_WORLD_BORDERS.shp | head -80" ] }, { "cell_type": "code", "execution_count": 32, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Default software rendering mode (use -h if accelerated video card installed).\n", "Loading tools from /usr/bin/openev/tools/Tool_Export.py\n", "Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py\n", "Loading tools from /usr/bin/openev/tools/Tool_DriverList.py\n" ] } ], "source": [ "%%bash\n", "# base on an attribute\n", "rm -f shp/TM_UGANDA_BORDERS-0.3.*\n", "ogr2ogr -overwrite -f \"ESRI Shapefile\" -where \"NAME = 'Uganda'\" geodata/shp/TM_UGANDA_BORDERS-0.3.shp geodata/shp/TM_WORLD_BORDERS.shp\n", "/usr/bin/openev/bin/openev geodata/shp/TM_UGANDA_BORDERS-0.3.shp" ] }, { "cell_type": "code", "execution_count": 33, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Default software rendering mode (use -h if accelerated video card installed).\n", "Loading tools from /usr/bin/openev/tools/Tool_Export.py\n", "Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py\n", "Loading tools from /usr/bin/openev/tools/Tool_DriverList.py\n" ] } ], "source": [ "%%bash\n", "# base on dimension of the polygons\n", "rm -f geodata/shp/TM_LARGE_BORDERS.*\n", "ogr2ogr -overwrite -f \"ESRI Shapefile\" -sql \"SELECT * FROM TM_WORLD_BORDERS WHERE OGR_GEOM_AREA > 100 \" geodata/shp/TM_LARGE_BORDERS.shp geodata/shp/TM_WORLD_BORDERS.shp\n", "/usr/bin/openev/bin/openev geodata/shp/TM_LARGE_BORDERS.shp" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Remove all the output" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [], "source": [ "%%bash\n", "rm -f geodata/shp/polygons.sqlite geodata/shp/point-buf_stat.csv geodata/shp/point_stat.csv geodata/vegetation/GPPcv08-11.tif LST/*_crop.tif geodata/vegetation/ETmean08-11_crop_trh.tif geodata/vegetation/ETmean08-11_crop_trh.txt vegetation/ETmean08-11_crop.txt geodata/vegetation/ETmosaic.vrt geodata/vegetation/ETmosaic.tif geodata/vegetation/stack_??.tif geodata/vegetation/stack.vrt geodata/vegetation/tiles.* geodata/vegetation/ETmean08-11_crop_proximity.tif geodata/vegetation/ETmean08-11_crop_buffer.tif geodata/dem/slope.tif geodata/dem/aspect.tif geodata/dem/tri.tif geodata/dem/tpi.tif geodata/dem/roughness.tif geodata/vegetation/ETmean08-11_01_trh?.tif geodata/LST/LST_MOYDmax_month1-msk.tif geodata/LST/LST_MOYDmax_monthStdev.tif geodata/LST/LST_MOYDmax_monthMean.tif geodata/LST/LST_MOYDmax_month_mean_stdev.tif geodata/LST/LST_MOYDmax_month.vrt geodata/LST/LST_MOYDmax_monthMean_aggreate10mean.tif geodata/LST/LST_MOYDmax_monthMean_circular11mean.tif geodata/temperature/reclass_ug_bio_3.tif geodata/temperature/reclass_ug_bio_3.txt geodata/shp/polygons_stat.csv geodata/shp/polygons_stat.* geodata/shp/TM_LARGE_BORDERS.shp.* geodata/shp/TM_UGANDA_BORDERS-0.3.* geodata/vegetation/ETmean08-11_crop.txt" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.8.10" } }, "nbformat": 4, "nbformat_minor": 4 }