Use GDAL/OGR for raster/vector operations - osgeo

Setting working directory for the OSGeoLive

GDAL 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/).

Open the bash terminal, migrate in the directory, and open the jupter-notebook.ipynb

cd $HOME/SE_data/exercise
jupyter-notebook 01_gdal_osgeo.ipynb

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.

Start to use GDAL commands

Explor the files

Change directory and list all the files with the extension .tif

[1]:
!ls $HOME/SE_data/exercise/geodata/*/*.tif
/home/user/SE_data/exercise/geodata/cloud/SA_intra.tif
/home/user/SE_data/exercise/geodata/cloud/SA_meanannual.tif
/home/user/SE_data/exercise/geodata/dem/GMTED2010.tif
/home/user/SE_data/exercise/geodata/dem/SA_elevation_mn_GMTED2010_mn.tif
/home/user/SE_data/exercise/geodata/landsat_ct/071W_41N_med_B1.tif
/home/user/SE_data/exercise/geodata/landsat_ct/071W_42N_med_B1.tif
/home/user/SE_data/exercise/geodata/landsat_ct/072W_40N_med_B1.tif
/home/user/SE_data/exercise/geodata/landsat_ct/072W_41N_med_B1.tif
/home/user/SE_data/exercise/geodata/landsat_ct/072W_42N_med_B1.tif
/home/user/SE_data/exercise/geodata/landsat_ct/073W_40N_med_B1.tif
/home/user/SE_data/exercise/geodata/landsat_ct/073W_41N_med_B1.tif
/home/user/SE_data/exercise/geodata/landsat_ct/073W_42N_med_B1.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month10.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month11.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month12.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month1.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month2.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month3.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month4.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month5.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month6.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month7.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month8.tif
/home/user/SE_data/exercise/geodata/LST/LST_MOYDmax_month9.tif
/home/user/SE_data/exercise/geodata/mask/msk_1km.tif
/home/user/SE_data/exercise/geodata/SDM/eBirdSampling_filtered_crop.tif
/home/user/SE_data/exercise/geodata/SDM/eBirdSampling_filtered.tif
/home/user/SE_data/exercise/geodata/temperature/ug_bio_3.tif
/home/user/SE_data/exercise/geodata/vegetation/ETmean08-11_01_msk.tif
/home/user/SE_data/exercise/geodata/vegetation/ETmean08-11.tif
/home/user/SE_data/exercise/geodata/vegetation/ETstdev08-11.tif
/home/user/SE_data/exercise/geodata/vegetation/GPPmean08-11.tif
/home/user/SE_data/exercise/geodata/vegetation/GPPstdev08-11.tif
/home/user/SE_data/exercise/geodata/vegetation/SA_tree_mn_percentage_GFC2013.tif

Retrive the characteristic of the one tif file

[2]:
!gdalinfo geodata/vegetation/ETmean08-11.tif
Driver: GTiff/GeoTIFF
Files: geodata/vegetation/ETmean08-11.tif
Size is 720, 600
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (29.000000000000000,4.000000004000000)
Pixel Size = (0.008333333340000,-0.008333333340000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  29.0000000,   4.0000000) ( 29d 0' 0.00"E,  4d 0' 0.00"N)
Lower Left  (  29.0000000,  -1.0000000) ( 29d 0' 0.00"E,  1d 0' 0.00"S)
Upper Right (  35.0000000,   4.0000000) ( 35d 0' 0.00"E,  4d 0' 0.00"N)
Lower Right (  35.0000000,  -1.0000000) ( 35d 0' 0.00"E,  1d 0' 0.00"S)
Center      (  32.0000000,   1.5000000) ( 32d 0' 0.00"E,  1d30' 0.00"N)
Band 1 Block=720x2 Type=Float32, ColorInterp=Gray
  NoData Value=-3.39999995214436425e+38

Visualize the image

[1]:
!  /usr/bin/openev/bin/openev geodata/vegetation/*.tif
Default software rendering mode (use -h if accelerated video card installed).
Loading tools from /usr/bin/openev/tools/Tool_DriverList.py
Loading tools from /usr/bin/openev/tools/Tool_Export.py
Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py
Reply to the following questions:
What is the pixel size?
How do you get min and max pixel values?

Understanding data type

Numerical System Decimal vs Binary

Ranges of GDAL data types

base^bit

Precision

Image Size

GDAL data type

Minimum

Maximum

Byte

0

255

2^8

39M

UInt16

0

65,535

2^16

78M

Int16, CInt16

-32,768

32,767

2^16

78M

UInt32

0

4,294,967,295

2^32

155M

Int32, CInt32

-2,147,483,648

2,147,483,647

2^32

155M

Float32, CFloat32

-3.4E38

3.4E38

2^38

Integer -16777216 16777216 (2^24)

155M

Float64, CFloat64

-1.79E308

1.79E308

2^308

309M

Understanding NoData Value

[4]:
!gdalinfo -mm geodata/vegetation/ETmean08-11.tif
Driver: GTiff/GeoTIFF
Files: geodata/vegetation/ETmean08-11.tif
Size is 720, 600
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (29.000000000000000,4.000000004000000)
Pixel Size = (0.008333333340000,-0.008333333340000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  29.0000000,   4.0000000) ( 29d 0' 0.00"E,  4d 0' 0.00"N)
Lower Left  (  29.0000000,  -1.0000000) ( 29d 0' 0.00"E,  1d 0' 0.00"S)
Upper Right (  35.0000000,   4.0000000) ( 35d 0' 0.00"E,  4d 0' 0.00"N)
Lower Right (  35.0000000,  -1.0000000) ( 35d 0' 0.00"E,  1d 0' 0.00"S)
Center      (  32.0000000,   1.5000000) ( 32d 0' 0.00"E,  1d30' 0.00"N)
Band 1 Block=720x2 Type=Float32, ColorInterp=Gray
    Computed Min/Max=0.547,7.492
  NoData Value=-3.39999995214436425e+38
[2]:
!gdal_edit.py -a_nodata -9999 geodata/vegetation/ETmean08-11.tif
[3]:
!gdalinfo -mm geodata/vegetation/ETmean08-11.tif
Driver: GTiff/GeoTIFF
Files: geodata/vegetation/ETmean08-11.tif
Size is 720, 600
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (29.000000000000000,4.000000004000000)
Pixel Size = (0.008333333340000,-0.008333333340000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  29.0000000,   4.0000000) ( 29d 0' 0.00"E,  4d 0' 0.00"N)
Lower Left  (  29.0000000,  -1.0000000) ( 29d 0' 0.00"E,  1d 0' 0.00"S)
Upper Right (  35.0000000,   4.0000000) ( 35d 0' 0.00"E,  4d 0' 0.00"N)
Lower Right (  35.0000000,  -1.0000000) ( 35d 0' 0.00"E,  1d 0' 0.00"S)
Center      (  32.0000000,   1.5000000) ( 32d 0' 0.00"E,  1d30' 0.00"N)
Band 1 Block=720x2 Type=Float32, ColorInterp=Gray
    Computed Min/Max=-339999995214436424907732413799364296704.000,7.492
  NoData Value=-9999
[5]:
!gdal_edit.py -a_nodata -339999995214436424907732413799364296704.000    geodata/vegetation/ETmean08-11.tif
[8]:
!gdalinfo -mm geodata/vegetation/ETmean08-11.tif
Driver: GTiff/GeoTIFF
Files: geodata/vegetation/ETmean08-11.tif
Size is 720, 600
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (29.000000000000000,4.000000004000000)
Pixel Size = (0.008333333340000,-0.008333333340000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=LZW
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  29.0000000,   4.0000000) ( 29d 0' 0.00"E,  4d 0' 0.00"N)
Lower Left  (  29.0000000,  -1.0000000) ( 29d 0' 0.00"E,  1d 0' 0.00"S)
Upper Right (  35.0000000,   4.0000000) ( 35d 0' 0.00"E,  4d 0' 0.00"N)
Lower Right (  35.0000000,  -1.0000000) ( 35d 0' 0.00"E,  1d 0' 0.00"S)
Center      (  32.0000000,   1.5000000) ( 32d 0' 0.00"E,  1d30' 0.00"N)
Band 1 Block=720x2 Type=Float32, ColorInterp=Gray
    Computed Min/Max=0.547,7.492
  NoData Value=-3.39999995214436425e+38

Calculate minimum and maximum values for all the images

[6]:
%%bash
for file in geodata/vegetation/*.tif ; do
gdalinfo  -mm  $file  | grep Computed
done
    Computed Min/Max=-10.000,5.000
    Computed Min/Max=0.547,7.492
    Computed Min/Max=0.209,3.332
    Computed Min/Max=0.000,9.738
    Computed Min/Max=0.000,4.220
    Computed Min/Max=0.000,10000.000
    Computed Min/Max=0.000,10000.000
    Computed Min/Max=0.000,10000.000

Create a Coefficient of variation

GPPmean08-11.tif temporal mean of Gross Primary Productivity
GPPstdev08-11.tif temporal standard deviation of Gross Primary Productivity
[7]:
%%bash
gdal_calc.py --type=Float32 --overwrite  -A geodata/vegetation/GPPstdev08-11.tif -B  geodata/vegetation/GPPmean08-11.tif \
--co=COMPRESS=DEFLATE --co=ZLEVEL=9   --calc="( A.astype(float) / ( B.astype(float) + 0.000000001 ) )" --outfile=geodata/vegetation/GPPcv08-11.tif
0 .. 10 .. 20 .. 30 .. 40 .. 50 .. 60 .. 70 .. 80 .. 90 .. 100 - Done
[8]:
! gdalinfo -mm geodata/vegetation/GPPcv08-11.tif
Driver: GTiff/GeoTIFF
Files: geodata/vegetation/GPPcv08-11.tif
Size is 720, 600
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (29.000000000000000,4.000000004000000)
Pixel Size = (0.008333333340000,-0.008333333340000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=DEFLATE
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  29.0000000,   4.0000000) ( 29d 0' 0.00"E,  4d 0' 0.00"N)
Lower Left  (  29.0000000,  -1.0000000) ( 29d 0' 0.00"E,  1d 0' 0.00"S)
Upper Right (  35.0000000,   4.0000000) ( 35d 0' 0.00"E,  4d 0' 0.00"N)
Lower Right (  35.0000000,  -1.0000000) ( 35d 0' 0.00"E,  1d 0' 0.00"S)
Center      (  32.0000000,   1.5000000) ( 32d 0' 0.00"E,  1d30' 0.00"N)
Band 1 Block=720x2 Type=Float32, ColorInterp=Gray
    Computed Min/Max=0.000,1.476
  NoData Value=3.40282346600000016e+38

Change pixel resolution

Looping trough the images

[1]:
%%bash
for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do
  filename=$(basename $file .tif )
  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
done
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month1.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month1.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month1.tif to destination geodata/LST/LST_MOYDmax_month1_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month2.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month2.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month2.tif to destination geodata/LST/LST_MOYDmax_month2_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month3.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month3.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month3.tif to destination geodata/LST/LST_MOYDmax_month3_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month4.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month4.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month4.tif to destination geodata/LST/LST_MOYDmax_month4_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month5.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month5.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month5.tif to destination geodata/LST/LST_MOYDmax_month5_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month6.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month6.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month6.tif to destination geodata/LST/LST_MOYDmax_month6_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month7.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month7.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month7.tif to destination geodata/LST/LST_MOYDmax_month7_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month8.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month8.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month8.tif to destination geodata/LST/LST_MOYDmax_month8_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month9.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month9.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month9.tif to destination geodata/LST/LST_MOYDmax_month9_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month10.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month10.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month10.tif to destination geodata/LST/LST_MOYDmax_month10_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month11.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month11.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month11.tif to destination geodata/LST/LST_MOYDmax_month11_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.
Creating output file that is 720P x 600L.
Processing geodata/LST/LST_MOYDmax_month12.tif [1/1] : 0Using internal nodata values (e.g. -9999) for image geodata/LST/LST_MOYDmax_month12.tif.
Copying nodata values from source geodata/LST/LST_MOYDmax_month12.tif to destination geodata/LST/LST_MOYDmax_month12_crop.tif.
...10...20...30...40...50...60...70...80...90...100 - done.

Get value at one pixel/line image-location

Looping trough the images

[13]:
%%bash
for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do
   gdallocationinfo  $file  1079 0
done
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 50.8218383789062
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 52.4714660644531
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 50.2210998535156
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 45.1343078613281
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 43.8952941894531
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 45.8667602539062
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 47.9639892578125
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 49.0910339355469
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 49.9898986816406
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 47.5579528808594
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 45.6031494140625
Report:
  Location: (1079P,0L)
  Band 1:
    Value: 47.5201110839844

Get value at lat/long image-location

Looping trough the images

[14]:
%%bash
for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do
   gdallocationinfo -valonly -geoloc $file  33.3 2
done
37.4764404296875
40.2584228515625
38.2513122558594
33.0213928222656
30.0732116699219
29.0838928222656
28.9105224609375
29.6008605957031
31.0305480957031
30.5758056640625
30.636962890625
33.5576171875
Get value at multiple lat/long image-location
[16]:
%%bash
# Create the lat long file
echo 32.5 2.5 > geodata/LST/x_y.txt
echo 31.1 2.1 >> geodata/LST/x_y.txt
# looping trough the images
for file in geodata/LST/LST_MOYDmax_month?.tif geodata/LST/LST_MOYDmax_month??.tif; do
   gdallocationinfo -valonly -geoloc $file < geodata/LST/x_y.txt
   echo ""
done
37.4022827148438
35.0345764160156

40.3694458007812
37.824951171875

38.5549926757812
36.6663208007812

32.7738952636719
32.6803283691406

29.0351257324219
30.0338134765625

29.025634765625
29.3355102539062

29.0183715820312
29.4237365722656

29.1660461425781
29.3502197265625

29.6674194335938
29.7001647949219

28.8177490234375
29.066650390625

29.2975463867188
29.015869140625

32.8586730957031
31.3289184570312

From Image to text and from txt to image

[18]:
%%bash
# from tif to 3 col txt XYZ
gdal_translate -of  XYZ   geodata/vegetation/ETmean08-11.tif geodata/vegetation/ETmean08-11_crop.txt
# masking
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
# back to tif
gdal_translate -co COMPRESS=DEFLATE -co ZLEVEL=9  -ot Byte geodata/vegetation/ETmean08-11_crop_trh.txt geodata/vegetation/ETmean08-11_crop_trh.tif
Input file size is 720, 600
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 720, 600
0...10...20...30...40...50...60...70...80...90...100 - done.

The use of VRT to stack images and tiling

[19]:
%%bash
# stck the tif
gdalbuildvrt -overwrite -separate geodata/vegetation/stack.vrt geodata/vegetation/ETmean08-11.tif geodata/vegetation/ETstdev08-11.tif
gdalinfo geodata/vegetation/stack.vrt
0...10...20...30...40...50...60...70...80...90...100 - done.
Driver: VRT/Virtual Raster
Files: geodata/vegetation/stack.vrt
       geodata/vegetation/ETmean08-11.tif
       geodata/vegetation/ETstdev08-11.tif
Size is 720, 600
Coordinate System is:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0],
    UNIT["degree",0.0174532925199433],
    AUTHORITY["EPSG","4326"]]
Origin = (29.000000000000000,4.000000004000000)
Pixel Size = (0.008333333340000,-0.008333333340000)
Corner Coordinates:
Upper Left  (  29.0000000,   4.0000000) ( 29d 0' 0.00"E,  4d 0' 0.00"N)
Lower Left  (  29.0000000,  -1.0000000) ( 29d 0' 0.00"E,  1d 0' 0.00"S)
Upper Right (  35.0000000,   4.0000000) ( 35d 0' 0.00"E,  4d 0' 0.00"N)
Lower Right (  35.0000000,  -1.0000000) ( 35d 0' 0.00"E,  1d 0' 0.00"S)
Center      (  32.0000000,   1.5000000) ( 32d 0' 0.00"E,  1d30' 0.00"N)
Band 1 Block=128x128 Type=Float32, ColorInterp=Undefined
  NoData Value=-3.39999995214436387e+38
Band 2 Block=128x128 Type=Float32, ColorInterp=Undefined
  NoData Value=-3.39999995214436387e+38
[20]:
%%bash
# split in tiles
gdal_translate -srcwin 0     0 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_UL.tif
gdal_translate -srcwin 0   300 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_UR.tif
gdal_translate -srcwin 360   0 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_LL.tif
gdal_translate -srcwin 360 300 360 300 geodata/vegetation/stack.vrt geodata/vegetation/stack_LR.tif

# recompose the image
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
gdal_translate -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/vegetation/ETmosaic.vrt geodata/vegetation/ETmosaic.tif
Input file size is 720, 600
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 720, 600
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 720, 600
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 720, 600
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
Input file size is 720, 600
0...10...20...30...40...50...60...70...80...90...100 - done.

Create shp border tiles

[21]:
%%bash
rm -f geodata/vegetation/tiles.*
gdaltindex  geodata/vegetation/tiles.shp  geodata/vegetation/stack_??.tif
Creating new index file...

Surface distance/buffer

[22]:
%%bash
# Continues distance surface
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
0...10...20...30...40...50...60...70...80...90...100 - done.
[23]:
%%bash
# Fix buffer surface
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
0...10...20...30...40...50...60...70...80...90...100 - done.

Topography variables

[25]:
%%bash
# calculate  slope
gdaldem slope -s 111120 -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/slope.tif

# calculate  apect
gdaldem aspect -zero_for_flat -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/aspect.tif

# calculate  Terrain Ruggedness Index TRI
gdaldem TRI -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/tri.tif

# calculate  Topographic Position Index TPI
gdaldem TPI -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/tpi.tif

# calculate  Roughness
gdaldem roughness -co COMPRESS=DEFLATE -co ZLEVEL=9 geodata/dem/GMTED2010.tif geodata/dem/roughness.tif
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.
0...10...20...30...40...50...60...70...80...90...100 - done.

Start to use OGR Commands

Select a polygon and create a new shape file

[3]:
%%bash
ogrinfo -al -geom=NO   geodata/shp/TM_WORLD_BORDERS.shp  | head -80
INFO: Open of `geodata/shp/TM_WORLD_BORDERS.shp'
      using driver `ESRI Shapefile' successful.

Layer name: TM_WORLD_BORDERS
Metadata:
  DBF_DATE_LAST_UPDATE=2021-02-09
Geometry: Polygon
Feature Count: 246
Extent: (-180.000000, -90.000000) - (180.000000, 83.623596)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
FIPS: String (2.0)
ISO2: String (2.0)
ISO3: String (3.0)
UN: Integer (3.0)
NAME: String (50.0)
AREA: Integer (7.0)
POP2005: Integer64 (10.0)
REGION: Integer (3.0)
SUBREGION: Integer (3.0)
LON: Real (8.3)
LAT: Real (7.3)
OGRFeature(TM_WORLD_BORDERS):0
  FIPS (String) = AC
  ISO2 (String) = AG
  ISO3 (String) = ATG
  UN (Integer) = 28
  NAME (String) = Antigua and Barbuda
  AREA (Integer) = 44
  POP2005 (Integer64) = 83039
  REGION (Integer) = 19
  SUBREGION (Integer) = 29
  LON (Real) = -61.783
  LAT (Real) = 17.078

OGRFeature(TM_WORLD_BORDERS):1
  FIPS (String) = AG
  ISO2 (String) = DZ
  ISO3 (String) = DZA
  UN (Integer) = 12
  NAME (String) = Algeria
  AREA (Integer) = 238174
  POP2005 (Integer64) = 32854159
  REGION (Integer) = 2
  SUBREGION (Integer) = 15
  LON (Real) = 2.632
  LAT (Real) = 28.163

OGRFeature(TM_WORLD_BORDERS):2
  FIPS (String) = AJ
  ISO2 (String) = AZ
  ISO3 (String) = AZE
  UN (Integer) = 31
  NAME (String) = Azerbaijan
  AREA (Integer) = 8260
  POP2005 (Integer64) = 8352021
  REGION (Integer) = 142
  SUBREGION (Integer) = 145
  LON (Real) = 47.395
  LAT (Real) = 40.430

OGRFeature(TM_WORLD_BORDERS):3
  FIPS (String) = AL
  ISO2 (String) = AL
  ISO3 (String) = ALB
  UN (Integer) = 8
  NAME (String) = Albania
  AREA (Integer) = 2740
  POP2005 (Integer64) = 3153731
  REGION (Integer) = 150
  SUBREGION (Integer) = 39
[4]:
%%bash
# base on an attribute
rm -f shp/TM_UGANDA_BORDERS-0.3.*
ogr2ogr  -overwrite  -f "ESRI Shapefile"  -where "NAME = 'Uganda'" geodata/shp/TM_UGANDA_BORDERS-0.3.shp geodata/shp/TM_WORLD_BORDERS.shp
/usr/bin/openev/bin/openev geodata/shp/TM_UGANDA_BORDERS-0.3.shp
Default software rendering mode (use -h if accelerated video card installed).
Loading tools from /usr/bin/openev/tools/Tool_DriverList.py
Loading tools from /usr/bin/openev/tools/Tool_Export.py
Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py
[5]:
%%bash
# base on dimension of the polygons
rm -f geodata/shp/TM_LARGE_BORDERS.*
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
/usr/bin/openev/bin/openev geodata/shp/TM_LARGE_BORDERS.shp
Default software rendering mode (use -h if accelerated video card installed).
Loading tools from /usr/bin/openev/tools/Tool_DriverList.py
Loading tools from /usr/bin/openev/tools/Tool_Export.py
Loading tools from /usr/bin/openev/tools/Tool_ShapesGrid.py

Remove all the output

[26]:
%%bash
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