An error occurred while loading the file. Please try again.
-
Le Roux Erwan authored
[exceeding] modify plot settings for the paper. put "imap" instead of "map" to avoid computing 2 times the models
e4b82417
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
from osgeo import gdal, ogr, osr
import numpy as np
import otb_numpy_proc
import otbApplication as otb
import os
def pixels_to_vector(selection_array, otb_pipeline, out_vector, out_key='out', var_dict=None):
geo_transform = otb_pipeline.GetImageMetaData(out_key)['GeoTransform']
epsg = int(otb_pipeline.GetImageMetaData(out_key)['ProjectionRef'].split('"')[-2])
srs = osr.SpatialReference()
srs.ImportFromEPSG(epsg)
drv = ogr.GetDriverByName('ESRI Shapefile')
ds = drv.CreateDataSource(out_vector)
ly = ds.CreateLayer(os.path.basename(os.path.splitext(out_vector)[0]), geom_type=ogr.wkbPoint, srs=srs)
ly.CreateField(ogr.FieldDefn('id', ogr.OFTInteger64))
ly.CreateField(ogr.FieldDefn('px', ogr.OFTInteger64))
ly.CreateField(ogr.FieldDefn('py', ogr.OFTInteger64))
if var_dict is not None:
for v in var_dict:
if 'float' in str(v['array'].dtype):
typ = ogr.OFTReal
else:
typ = ogr.OFTInteger
if v['array'].ndim == 2:
for t in range(v['array'].shape[1]):
ly.CreateField(ogr.FieldDefn(v['prefix']+'_%d' % t, typ))
else:
ly.CreateField(ogr.FieldDefn(v['prefix'], typ))
i = 0
pix_pnts = np.where(selection_array)
for p in zip(pix_pnts[0]+0.5,pix_pnts[1]+0.5):
X = geo_transform[0] + p[1] * geo_transform[1]
Y = geo_transform[3] + p[0] * geo_transform[5]
gp = ogr.Geometry(ogr.wkbPoint)
gp.AddPoint(X,Y)
gf = ogr.Feature(ly.GetLayerDefn())
gf.SetGeometry(gp)
gf.SetField('id', i)
gf.SetField('px', p[1])
gf.SetField('py', p[0])
if var_dict is not None:
for v in var_dict:
if v['array'].ndim == 2:
for t in range(v['array'].shape[1]):
gf.SetField(v['prefix']+'_%d' % t, float(v['array'][i,t]))
else:
gf.SetField(v['prefix'], float(v['array'][i]))
ly.CreateFeature(gf)
i += 1
ds = None
if __name__ == '__main__':
p = otb_numpy_proc.to_otb_pipeline('./pippo.tif')
p.Execute()
arr = np.random.random((100, 100)) < 0.1
drv = ogr.GetDriverByName('ESRI Shapefile')
drv.DeleteDataSource('./pippuzzo.shp')
pixels_to_vector(arr, p, 'pippuzzo.shp')