Geomatics/gdal-pointInTiff.py at master · yoba7/Geomatics · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Sun Feb  7 00:09:07 2021

@author: yoba
"""

import numpy as np
from osgeo import gdal
from numpy.linalg import inv as inverse
import pandas as pd


'''
  Parameters:

'''

tifZipFile='../ref/belgium_clc.zip'
tifFile='belgium_clc.tif'
fileContainingXYCoordinates='../input/sample-from-best-3035.zip'
fileContainingXYCoordinatesMergedWithTiff='../output/sample-from-best-3035-with-clc.zip'
compression_opts = dict(method='zip', archive_name='sample-from-best-3035-with-clc.csv')

tifZipFile='../ref/belgium_dem_v11.zip'
tifFile='belgium_dem_v11.tif'
fileContainingXYCoordinates='../input/sample-from-best-3035.zip'
fileContainingXYCoordinatesMergedWithTiff='../output/sample-from-best-3035-with-dem.zip'
compression_opts = dict(method='zip', archive_name='sample-from-best-3035-with-dem.csv')

'''
   Use gdal's virtual file system to access a tif inside a zip.

   See: https://gdal.org/user/virtual_file_systems.html
'''


myRaster = gdal.Open(f'/vsizip/{tifZipFile}/{tifFile}')
#myRaster = gdal.Open(f'{tifFile}')


'''
 "Connect" to first band and tranform it into a matrix.

 There is only one band in our Corine Land Cover file, but other
 geotiffs might have multiple bands.

 The object returned by ReadAsAssay() is a numpy array:

 >> type(myBand)
 >> numpy.ndarray

 The numpy array can be very large. You can get the size of the array
 using myBand.shape. In this example, the size of the matry is:

 >> myBand.shape
 >> (2752, 3247)


'''

myBand=myRaster.GetRasterBand(1).ReadAsArray()

'''
 Extract metadata from raster file.

 See: https://gdal.org/user/raster_data_model.html#affine-geotransform

 The metadata makes it possible to create a matrix to
 transform (pixel,line) coordinates into (x,y) georeferenced space.

 In our example, the pixel2xy matrix will transform (pixel,line) coordinates
 of belgium_clc.tif into (x,y) coordinates in the EPSG:3035 projection system.

'''

M=myRaster.GetGeoTransform()

pixel2xy=np.array([[M[1], M[4]],
                   [M[2], M[5]]])

translation=np.array([M[0],M[3]])


'''
 Compute "xy2pixel", ie, the inverse affine transformation matrix.

 This "xy2pixel" matrix can be used to convert (x,y) coordinates
 into (pixel,line) coordinates.

'''

xy2pixel=inverse(pixel2xy)

'''
 Import a file containing (x,y) coordinates into a Pandas dataframe.
'''

adresses=pd.read_csv(f'{fileContainingXYCoordinates}',sep=',')

'''
 Convert into numpy array (matrix) and apply transformation
'''

xy=np.array(adresses[['x','y']])

pl = (xy-translation) @ xy2pixel

pl=pd.DataFrame(pl,columns=['p','l'])

lmax,pmax=myBand.shape

'''
 Extract information from matrix
'''

def extractInfoFor(l,p):
    if l>=0 and l<lmax and p>=0 and p<pmax:
        return myBand[int(l),int(p)]
    else:
        return -9999

adresses['value']=pl.apply(lambda row: extractInfoFor(row.l,row.p),axis=1)

'''
 Export results
'''

adresses.to_csv(f'{fileContainingXYCoordinatesMergedWithTiff}',
                sep='|',
                index=False,
                compression=compression_opts)