-
Notifications
You must be signed in to change notification settings - Fork 7
/
index_setsm.py
executable file
·1360 lines (1164 loc) · 66.3 KB
/
index_setsm.py
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
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
import argparse
import configparser
import datetime
import json
import logging
import os
import pickle
import re
import sys
from osgeo import gdal, osr, ogr
from lib import utils, dem, walk
from lib import VERSION, SHORT_VERSION
logger = utils.get_logger()
utils.setup_gdal_error_handler()
gdal.UseExceptions()
# Script paths and execution
SCRIPT_FILE = os.path.abspath(os.path.realpath(__file__))
SCRIPT_FNAME = os.path.basename(SCRIPT_FILE)
SCRIPT_NAME, SCRIPT_EXT = os.path.splitext(SCRIPT_FNAME)
SCRIPT_DIR = os.path.dirname(SCRIPT_FILE)
FORMAT_OPTIONS = {
'SHP':'ESRI Shapefile',
'GDB':'ESRI Geodatabase',
'PG':'PostgreSQL Database (PG:<config.ini section with connection info>:<layer name>}',
}
FORMAT_HELP = ['{}:{},'.format(k,v) for k, v in FORMAT_OPTIONS.items()]
PROJECTS = {
'arcticdem': 'ArcticDEM',
'rema': 'REMA',
'earthdem': 'EarthDEM',
}
mask_strip_suffixes = (
'_dem_water-masked.tif',
'_dem_cloud-masked.tif',
'_dem_cloud-water-masked.tif',
'_dem_masked.tif'
)
MODES = {
## mode : (class, suffix, groupid_fld, field_def)
'scene': (dem.SetsmScene, '_meta.txt', 'stripdemid',
utils.SCENE_ATTRIBUTE_DEFINITIONS, utils.SCENE_ATTRIBUTE_DEFINITIONS_REGISTRATION),
'strip': (dem.SetsmDem, '_dem.tif', 'stripdirname',
utils.DEM_ATTRIBUTE_DEFINITIONS, utils.DEM_ATTRIBUTE_DEFINITIONS_REGISTRATION),
'tile': (dem.SetsmTile, '_dem.tif', 'supertile_id',
utils.TILE_DEM_ATTRIBUTE_DEFINITIONS, utils.TILE_DEM_ATTRIBUTE_DEFINITIONS_REGISTRATION),
}
id_flds = ['SCENEDEMID', 'STRIPDEMID', 'DEM_ID', 'TILE', 'LOCATION', 'INDEX_DATE', 'IS_DSP']
recordid_map = {
'scene': '{SCENEDEMID}|{STRIPDEMID}|{IS_DSP}|{LOCATION}|{INDEX_DATE}',
'strip': '{DEM_ID}|{STRIPDEMID}|{LOCATION}|{INDEX_DATE}',
'strip_release': '{DEM_ID}|{STRIPDEMID}|{FILEURL}|{CR_DATE}',
'tile': '{DEM_ID}|{TILE}|{LOCATION}|{INDEX_DATE}',
'tile_release': '{DEM_ID}|{TILE}|{FILEURL}|{CR_DATE}',
}
BP_PATH_PREFIX = 'https://blackpearl-data2.pgc.umn.edu'
PGC_PATH_PREFIX = '/mnt/pgc/data/elev/dem/setsm'
BW_PATH_PREFIX = '/scratch/sciteam/GS_bazu/elev/dem/setsm'
CSS_PATH_PREFIX = '/css/nga-dems/setsm'
custom_path_prefixes = {
'BP': BP_PATH_PREFIX,
'PGC': PGC_PATH_PREFIX,
'BW': BW_PATH_PREFIX,
'CSS': CSS_PATH_PREFIX
}
DSP_OPTIONS = {
'dsp': 'record using current downsample product DEM res',
'orig': 'record using original pre-DSP DEM res',
'both': 'write a record for each'
}
dem_type_folder_lookup = {
'strip': 'strips',
'tile': 'mosaics',
}
DEFAULT_DSP_OPTION = 'dsp'
# handle unicode in Python 3
try:
unicode('')
except NameError:
unicode = str
class RawTextArgumentDefaultsHelpFormatter(argparse.ArgumentDefaultsHelpFormatter, argparse.RawTextHelpFormatter): pass
def main():
#### Set Up Arguments
parser = argparse.ArgumentParser(
formatter_class=RawTextArgumentDefaultsHelpFormatter,
description="build setsm DEM index"
)
#### Positional Arguments
parser.add_argument('src', help="source directory or image")
parser.add_argument('dst', help="destination index dataset (use PG:<config.ini section name>:<layer name> for a postgresql DB")
#### Optional Arguments
parser.add_argument('--mode', choices=MODES.keys(), default='scene',
help="type of items to index {} default=scene".format(MODES.keys()))
parser.add_argument('--config', default=os.path.join(SCRIPT_DIR, 'config.ini'),
help="config file (default is config.ini in script dir, or fallback to ~/.pg_service.conf)")
parser.add_argument('--epsg', type=int, default=4326,
help="egsg code for output index projection (default wgs85 geographic epsg:4326)")
parser.add_argument('--dsp-record-mode', choices=DSP_OPTIONS.keys(), default=DEFAULT_DSP_OPTION,
help='resolution mode for downsampled product (dsp) record (mode=scene only):\n{}'.format(
'\n'.join([k+': '+v for k,v in DSP_OPTIONS.items()])
))
parser.add_argument('--status', help='custom value for status field')
parser.add_argument('--status-dsp-record-mode-orig', help='custom value for status field when dsp-record-mode is set to "orig"')
# DEPRECATED
# parser.add_argument('--include-registration', action='store_true', default=False,
# help='include registration info if present (mode=strip and tile only)')
parser.add_argument('--use-release-fields', action='store_true', default=False,
help="use field definitions for tile release indices (mode=tile only)")
parser.add_argument('--long-fieldnames', action='store_true', default=False,
help="use long format (>10 chars) version of fieldnames")
parser.add_argument('--lowercase-fieldnames', action='store_true', default=False,
help="make fieldnames lowercase when writing to new destination index")
parser.add_argument('--search-masked', action='store_true', default=False,
help='search for masked and unmasked DEMs (mode=strip only)')
parser.add_argument('--read-json', action='store_true', default=False,
help='search for json files instead of images to populate the index')
parser.add_argument('--write-json', action='store_true', default=False,
help='write results to json files in dst folder')
parser.add_argument('--maxdepth', type=float, default=float('inf'),
help='maximum depth into source directory to be searched')
parser.add_argument('--log', help="directory for log output (debug messages written here)")
parser.add_argument('--overwrite', action='store_true', default=False,
help="overwrite existing index")
parser.add_argument('--append', action='store_true', default=False,
help="append records to existing index")
parser.add_argument('--check', action='store_true', default=False,
help='verify new records exist in target index (not compatible with --write-json or --dryrun)')
parser.add_argument('--skip-region-lookup', action='store_true', default=False,
help="skip region lookup on danco")
parser.add_argument('--skip-records-missing-dsp-original-info', action='store_true', default=False,
help="skip adding records where the file info on the source DEM for a dsp product is missing"
" (valid only if --dsp-record-mode is orig or both)")
parser.add_argument("--write-pickle", help="store region lookup in a pickle file. skipped if --write-json is used")
parser.add_argument("--read-pickle", help='read region lookup from a pickle file. skipped if --write-json is used')
parser.add_argument("--custom-paths", choices=custom_path_prefixes.keys(), help='Use custom path schema')
parser.add_argument('--project', choices=utils.PROJECTS.keys(), help='project name (required when writing tiles)')
parser.add_argument('--debug', action='store_true', default=False, help='print DEBUG level logger messages to terminal')
parser.add_argument('--dryrun', action='store_true', default=False, help='run script without inserting records')
parser.add_argument('--np', action='store_true', default=False, help='do not print progress bar')
parser.add_argument('--version', action='version', version=f"Current version: {SHORT_VERSION}",
help='print version and exit')
parser.add_argument('--release-fileurl', type=str, default="https://data.pgc.umn.edu/elev/dem/setsm/<project>/<type>/<version>/<resolution>/<group>/<dem_id>.tar.gz",
help="template for release field 'fileurl' (--use-release-fields only)")
parser.add_argument('--release-s3url', type=str, default="https://polargeospatialcenter.github.io/stac-browser/#/external/pgc-opendata-dems.s3.us-west-2.amazonaws.com/<project>/<type>/<version>/<resolution>/<group>/<dem_id>.json",
help="template for release field 's3url' (--use-release-fields only)")
#### Parse Arguments
args = parser.parse_args()
if args.debug:
utils.logger_streamhandler_debug()
#### Verify Arguments
if not os.path.isdir(args.src) and not os.path.isfile(args.src):
parser.error("Source directory or file does not exist: %s" %args.src)
src = args.src
dst = args.dst
if args.overwrite and args.append:
parser.error('--append and --overwrite are mutually exclusive')
if args.write_json and args.append:
parser.error('--append cannot be used with the --write-json option')
if (args.write_json or args.read_json) and args.search_masked:
parser.error('--search-masked cannot be used with the --write-json or --read-json options')
if args.mode != 'strip' and args.search_masked:
parser.error('--search-masked applies only to mode=strip')
if args.write_json and args.check:
parser.error('--check cannot be used with the --write-json option')
## Check project
if args.mode == 'tile' and not args.project:
parser.error("--project option is required if when mode=tile")
if args.mode == 'strip' and args.use_release_fields and not args.project:
parser.error("--project option is required when mode=strip using --use-release-fields")
if args.mode == 'scene' and args.use_release_fields:
parser.error("--use-release-fields option is not applicable to mode=scene")
## Todo add Bp region lookup via API instead of Danco?
if args.skip_region_lookup and (args.custom_paths == 'PGC' or args.custom_paths == 'BP'):
parser.error('--skip-region-lookup is not compatible with --custom-paths = PGC or BP')
if args.write_pickle:
if not os.path.isdir(os.path.dirname(args.write_pickle)):
parser.error("Pickle file must be in an existing directory")
if args.read_pickle:
if not os.path.isfile(args.read_pickle):
parser.error("Pickle file must be an existing file")
if args.status and args.custom_paths == 'BP':
parser.error("--custom_paths BP sets status field to 'tape' and cannot be used with --status. For dsp-record-mode=orig custom status, use --status-dsp-record-mode-orig")
path_prefix = custom_path_prefixes[args.custom_paths] if args.custom_paths else None
if args.log:
if os.path.isdir(args.log):
tm = datetime.datetime.now()
logfile = os.path.join(args.log,"index_setsm_{}.log".format(tm.strftime("%Y%m%d%H%M%S")))
else:
parser.error('log folder does not exist: {}'.format(args.log))
lfh = logging.FileHandler(logfile)
lfh.setLevel(logging.DEBUG)
formatter = logging.Formatter('%(asctime)s %(levelname)s- %(message)s','%m-%d-%Y %H:%M:%S')
lfh.setFormatter(formatter)
logger.addHandler(lfh)
## Check path if --write-json is invoked
if args.write_json:
if not os.path.isdir(dst):
parser.error("Destination must be an existing directory with --write-json option")
ogr_driver_str = None
if args.epsg:
logger.warning('--epsg and --dsp-original-res will be ignored with the --write-json option')
logger.info("Current repo version: %s", VERSION)
rc = 0
if args.write_json:
logger.info("Forcing indexer to use absolute paths for writing JSONs")
src = os.path.abspath(args.src)
## If not writing to JSON, get OGR driver, ds name, and layer name
else:
try:
ogr_driver_list, dst_dsp, dst_lyr = utils.get_source_names2(dst)
except RuntimeError as e:
parser.error(e)
ogrDriver = None
for ogr_driver_str in ogr_driver_list:
ogrDriver = ogr.GetDriverByName(ogr_driver_str)
if ogrDriver is not None:
break
if ogrDriver is None:
parser.error("Driver(s) not available: {}".format(', '.join(ogr_driver_list)))
else:
logger.info("Driver selected: {}".format(ogr_driver_str))
#### Get Config file contents
config = configparser.ConfigParser()
config.read(args.config)
pg_config_file = os.path.expanduser("~/.pg_service.conf")
pg_config = configparser.ConfigParser()
pg_config.read(pg_config_file)
#### Get output DB connection if specified
if ogr_driver_str in ("PostgreSQL"):
section = dst_dsp
if section in config.sections():
conn_info = {
'host':config.get(section,'host'),
'port':config.getint(section,'port'),
'name':config.get(section,'name'),
'schema':config.get(section,'schema'),
'user':config.get(section,'user'),
'pw':config.get(section,'pw'),
}
dst_ds = "PG:host={host} port={port} dbname={name} user={user} password={pw} active_schema={schema}".format(**conn_info)
conn_str_redacted = re.sub(r"password=\S+", "password=PASS", dst_ds)
logger.info(f"Derived dst dataset PG connection string from {args.config}: '{conn_str_redacted}'")
elif section in pg_config.sections():
dst_ds = f"PG:service={section}"
logger.info(f"Derived dst dataset PG connection string from {pg_config_file}: '{dst_ds}'")
else:
logger.error(f"--config file or ~/.pg_service.conf must contain credentials for service name '{section}'")
rc = -1
#### Set dataset path is SHP or GDB
elif ogr_driver_str in ("ESRI Shapefile", "FileGDB", "OpenFileGDB", "GPKG"):
dst_ds = dst_dsp
if not os.path.isdir(os.path.dirname(dst_dsp)):
parser.error('--dst must be within an existing directory')
else:
logger.error("Format {} is not supported".format(ogr_driver_str))
rc = -1
## Get pairname-region dict
if args.skip_region_lookup or args.mode == 'tile':
pairs = {}
else:
if args.read_pickle:
logger.info("Fetching region lookup from pickle file")
pairs = pickle.load(open(args.read_pickle, "rb"))
else:
#### Get Danco connection if available
section_depr = 'danco'
section = 'pgc_danco_footprint'
conn_str = None
if section not in config.sections() and section_depr in config.sections():
logger.warning(f"Config section name '{section_depr}' is deprecated and should be changed to '{section}'")
section = section_depr
if section in config.sections():
danco_conn_info = {
'host':config.get(section,'host'),
'port':config.getint(section,'port'),
'name':config.get(section,'name'),
'schema':config.get(section,'schema'),
'user':config.get(section,'user'),
'pw':config.get(section,'pw'),
}
conn_str = "PG:host={host} port={port} dbname={name} user={user} password={pw} active_schema={schema}".format(**danco_conn_info)
conn_str_redacted = re.sub(r"password=\S+", "password=PASS", conn_str)
logger.info(f"Derived Danco connection string from {args.config}: '{conn_str_redacted}'")
elif section in pg_config.sections():
conn_str = f"PG:service={section} active_schema=public"
logger.info(f"Derived Danco connection string from {pg_config_file}: '{conn_str}'")
if conn_str:
logger.info("Fetching region lookup from Danco")
pairs = get_pair_region_dict(conn_str)
else:
logger.warning(f"--config file or ~/.pg_service.conf do not contain credentials for service name '{section}'. Region cannot be determined.")
pairs = {}
if len(pairs) == 0:
logger.warning("Cannot get region-pair lookup")
if args.custom_paths == 'PGC' or args.custom_paths == 'BP':
logger.error("Region-pair lookup required for --custom_paths PGC or BP option")
sys.exit()
## Save pickle if selected
if args.write_pickle:
logger.info("Pickling region lookup")
pickle.dump(pairs, open(args.write_pickle, "wb"))
#### Test epsg
try:
spatial_ref = utils.SpatialRef(args.epsg)
except RuntimeError as e:
parser.error(e)
#### Test if dst table exists
if ogr_driver_str == 'ESRI Shapefile':
if os.path.isfile(dst_ds):
if args.overwrite:
logger.info("Removing old index... %s" %os.path.basename(dst_ds))
if not args.dryrun:
ogrDriver.DeleteDataSource(dst_ds)
elif not args.append:
logger.error("Dst shapefile exists. Use the --overwrite or --append options.")
rc = -1
elif ogr_driver_str in ('FileGDB', 'OpenFileGDB', 'GPKG'):
if os.path.isdir(dst_ds):
ds = ogrDriver.Open(dst_ds,1)
if ds:
for i in range(ds.GetLayerCount()):
lyr = ds.GetLayer(i)
if lyr.GetName() == dst_lyr:
if args.overwrite:
logger.info("Removing old index layer: {}".format(dst_lyr))
del lyr
ds.DeleteLayer(i)
break
elif not args.append:
logger.error("Dst GDB layer exists. Use the --overwrite or --append options.")
rc = -1
ds = None
## Postgres check - do not overwrite
elif ogr_driver_str == 'PostgreSQL':
ds = ogrDriver.Open(dst_ds,1)
if ds:
for i in range(ds.GetLayerCount()):
lyr = ds.GetLayer(i)
if lyr.GetName() == dst_lyr:
if args.overwrite:
logger.info("Removing old index layer: {}".format(dst_lyr))
del lyr
ds.DeleteLayer(i)
break
elif not args.append:
logger.error("Dst DB layer exists. Use the --overwrite or --append options.")
rc = -1
ds = None
else:
logger.error("Format {} not handled in dst table existence check".format(ogr_driver_str))
rc = -1
if rc == 0:
#### ID records
dem_class, suffix, groupid_fld, fld_defs_base, reg_fld_defs = MODES[args.mode]
if args.mode == 'tile' and args.use_release_fields:
fld_defs_base = utils.TILE_DEM_ATTRIBUTE_DEFINITIONS_RELEASE
if args.mode == 'strip' and args.use_release_fields:
fld_defs_base = utils.DEM_ATTRIBUTE_DEFINITIONS_RELEASE
if args.mode == 'strip' and args.search_masked:
suffix = mask_strip_suffixes + tuple([suffix])
# fld_defs = fld_defs_base + reg_fld_defs if args.include_registration else fld_defs_base - DEPRECATED
fld_defs = fld_defs_base
src_fps = []
records = []
logger.info('Source: {}'.format(src))
logger.info('Identifying DEMs')
if os.path.isfile(src):
logger.info(src)
src_fps.append(src)
else:
for root, dirs, files in walk.walk(src, maxdepth=args.maxdepth):
for f in files:
if (f.endswith('.json') and args.read_json) or (f.endswith(suffix) and not args.read_json):
logger.debug(os.path.join(root,f))
src_fps.append(os.path.join(root,f))
total = len(src_fps)
i=0
for src_fp in src_fps:
i+=1
if not args.np:
utils.progress(i, total, "DEMs identified")
if args.read_json:
temp_records = read_json(os.path.join(src_fp),args.mode)
records.extend(temp_records)
else:
record = None
try:
record = dem_class(src_fp)
record.get_dem_info()
except Exception as e:
logger.error(e)
if record is not None and hasattr(record, 'srcfp'):
logger.error("Error encountered on DEM record: {}".format(record.srcfp))
else:
## Check if DEM is a DSP DEM, dsp-record mode includes 'orig', and the original DEM data is unavailable
if args.mode == 'scene' and record.is_dsp and not os.path.isfile(record.dspinfo) \
and args.dsp_record_mode in ['orig', 'both']:
logger.error("DEM {} has no Dsp downsample info file: {}, skipping".format(record.id,record.dspinfo))
else:
records.append(record)
if not args.np:
print('')
total = len(records)
if total == 0:
logger.info("No valid records found")
else:
logger.info("{} records found".format(total))
## Group into strips or tiles for json writing
groups = {}
for record in records:
groupid = getattr(record, groupid_fld)
if groupid in groups:
groups[groupid].append(record)
else:
groups[groupid] = [record]
#### Write index
if args.write_json:
write_to_json(dst, groups, total, args)
elif not args.dryrun:
rc = write_to_ogr_dataset(ogr_driver_str, ogrDriver, dst_ds, dst_lyr, groups,
pairs, total, path_prefix, fld_defs, args)
else:
logger.info("Exiting dryrun")
sys.exit(rc)
def write_to_ogr_dataset(ogr_driver_str, ogrDriver, dst_ds, dst_lyr, groups, pairs, total, path_prefix, fld_defs, args):
ds = None
rc = 0
## Create dataset if it does not exist
if ogr_driver_str == 'ESRI Shapefile':
max_fld_width = 254
if os.path.isfile(dst_ds):
ds = ogrDriver.Open(dst_ds,1)
else:
ds = ogrDriver.CreateDataSource(dst_ds)
elif ogr_driver_str in ['FileGDB', 'OpenFileGDB', 'GPKG']:
max_fld_width = 1024
if os.path.isdir(dst_ds):
ds = ogrDriver.Open(dst_ds,1)
else:
ds = ogrDriver.CreateDataSource(dst_ds)
elif ogr_driver_str == 'PostgreSQL':
max_fld_width = 1024
# DB must already exist
ds = ogrDriver.Open(dst_ds,1)
else:
logger.error("Format {} is not supported".format(ogr_driver_str))
if args.status:
status = args.status
elif args.custom_paths == 'BP':
status = 'tape'
else:
status = 'online'
dsp_orig_status = args.status_dsp_record_mode_orig if args.status_dsp_record_mode_orig else status
fld_def_location_fwidth_gdb = None
for f in fld_defs:
if f.fname.upper() == 'LOCATION':
fld_def_location_fwidth_gdb = min(f.fwidth, 1024)
break
fld_def_short_to_long_dict = {
field_def.fname: (field_def.fname_long if field_def.fname_long else field_def.fname) for field_def in fld_defs
}
if ds is None:
logger.info("Cannot open dataset: {}".format(dst_ds))
rc = -1
else:
## Create table if it does not exist
layer = ds.GetLayerByName(dst_lyr)
fld_list = [f.fname for f in fld_defs]
tgt_srs = utils.osr_srs_preserve_axis_order(osr.SpatialReference())
tgt_srs.ImportFromEPSG(args.epsg)
if not layer:
logger.info("Creating table...")
# FileGDB will throw a warning when inserting datetimes without this
if ogr_driver_str in ['FileGDB', 'OpenFileGDB']:
co = ['TIME_IN_UTC=NO']
else:
co = []
layer = ds.CreateLayer(dst_lyr, tgt_srs, ogr.wkbMultiPolygon, options=co)
if layer:
for field_def in fld_defs:
fname = fld_def_short_to_long_dict[field_def.fname] if args.long_fieldnames else field_def.fname
if args.lowercase_fieldnames:
fname = fname.lower()
fstype = None
if field_def.ftype == ogr.OFTDateTime and ogr_driver_str in ['ESRI Shapefile']:
ftype = ogr.OFTString
fwidth = 28
elif field_def.ftype == ogr.OFSTBoolean:
ftype = ogr.OFTInteger
fstype = field_def.ftype
fwidth = field_def.fwidth
else:
ftype = field_def.ftype
fwidth = field_def.fwidth
field = ogr.FieldDefn(fname, ftype)
if fstype:
field.SetSubType(fstype)
field.SetWidth(min(max_fld_width, fwidth))
field.SetPrecision(field_def.fprecision)
layer.CreateField(field)
# When creating a new dataset/layer schema with GDAL, something about the schema
# is kept in cache tied to the dataset connection instance.
# If after writing to the new layer, the layer records are read back in with GDAL
# using the same connection instance (with the --check option for example),
# GDAL may fail to read certain field data types properly (such as boolean fields).
# This is likely a bug in GDAL!
# To avoid this, manually close and reopen the connection.
ds = None
layer = None
ds = ogrDriver.Open(dst_ds, 1)
layer = ds.GetLayerByName(dst_lyr)
## Append Records
if layer:
# Get field widths
lyr_def = layer.GetLayerDefn()
fwidths = {lyr_def.GetFieldDefn(i).GetName():
(lyr_def.GetFieldDefn(i).GetWidth(), lyr_def.GetFieldDefn(i).GetType())
for i in range(lyr_def.GetFieldCount())
}
fnameupper_fnamelayer_dict = {k.upper(): k for k, v in fwidths.items()}
logger.info("Appending records...")
#### loop through records and add features
i=0
recordids = []
invalid_record_cnt = 0
duplicate_record_cnt = 0
dsp_modes = ['orig','dsp'] if args.dsp_record_mode == 'both' else [args.dsp_record_mode]
for groupid in groups:
for record in groups[groupid]:
for dsp_mode in dsp_modes:
region = None
bp_region = None
# skip writing a second "orig" record if the DEM is not a DSP DEM sene
if args.mode == 'scene':
if not record.is_dsp and dsp_mode == 'orig':
continue
i+=1
if not args.np:
utils.progress(i, total * len(dsp_modes), "features written")
feat = ogr.Feature(layer.GetLayerDefn())
valid_record = True
## Set attributes
## Fields for scene DEM
if args.mode == 'scene':
logger.debug(f"Processing scene: {record.sceneid} - mode {dsp_mode}")
attrib_map = {
'SCENEDEMID': record.dsp_sceneid if (dsp_mode == 'orig') else record.sceneid,
'STRIPDEMID': record.dsp_stripdemid if (dsp_mode == 'orig') else record.stripdemid,
'STATUS': dsp_orig_status if (dsp_mode == 'orig') else status,
'PAIRNAME': record.pairname,
'SENSOR1': record.sensor1,
'SENSOR2': record.sensor2,
'ACQDATE1': record.acqdate1.strftime('%Y-%m-%dT%H:%M:%SZ'),
'ACQDATE2': record.acqdate2.strftime('%Y-%m-%dT%H:%M:%SZ'),
'CATALOGID1': record.catid1,
'CATALOGID2': record.catid2,
'SCENE1': record.scene1,
'SCENE2': record.scene2,
'GEN_TIME1': record.gentime1.strftime('%Y-%m-%dT%H:%M:%SZ') if record.gentime1 else None,
'GEN_TIME2': record.gentime2.strftime('%Y-%m-%dT%H:%M:%SZ') if record.gentime2 else None,
'HAS_LSF': record.has_lsf,
'HAS_NONLSF': record.has_nonlsf,
'IS_XTRACK': record.is_xtrack,
'IS_DSP': False if dsp_mode == 'orig' else record.is_dsp,
'ALGM_VER': record.algm_version,
'PROD_VER': record.prod_version,
'PROJ4': record.proj4,
'EPSG': record.epsg,
}
attr_pfx = 'dsp_' if dsp_mode == 'orig' else ''
for k in record.filesz_attrib_map:
attrib_map[k.upper()] = getattr(record,'{}{}'.format(attr_pfx,k))
# TODO revisit after all incorrect 50cminfo.txt files are ingested
# Overwrite original res dsp filesz values will Null
if dsp_mode == 'orig':
for k in record.filesz_attrib_map:
attrib_map[k.upper()] = None
# Test if filesz attr is valid for dsp original res records
if dsp_mode == 'orig':
if attrib_map['FILESZ_DEM'] is None:
if not args.skip_records_missing_dsp_original_info:
logger.debug(
"Original res filesz_dem is empty for {}. Record will still be written".format(
record.sceneid))
else:
logger.error(
"Original res filesz_dem is empty for {}. Record skipped".format(record.sceneid))
valid_record = False
elif attrib_map['FILESZ_DEM'] == 0:
logger.warning(
"Original res filesz_dem is 0 for {}. Record will still be written".format(record.sceneid))
# Test if filesz attr is valid for normal records
elif not attrib_map['FILESZ_DEM'] and not attrib_map['FILESZ_LSF']:
logger.warning(
"DEM and LSF DEM file size is zero or null for {}. Record skipped".format(record.sceneid))
valid_record = False
# Set region
try:
pair_items = pairs[record.pairname]
except KeyError as e:
region = None
bp_region = None
else:
if isinstance(pair_items, str):
region = pair_items
elif isinstance(pair_items, tuple):
region, bp_region = pair_items
else:
logger.error("Pairname-region lookup value cannot be parsed for pairname {}: {}".format(
record.pairname, pair_items))
attrib_map['REGION'] = region
if path_prefix:
res_dir = record.res_str + '_dsp' if record.is_dsp else record.res_str
if args.custom_paths == 'BP':
# https://blackpearl-data2.pgc.umn.edu/dem-scenes-2m-arceua/2m/WV02/2015/05/
# WV02_20150506_1030010041510B00_1030010043050B00_50cm_v040002.tar
if not region:
logger.error("Pairname not found in region lookup {}, cannot build custom path".format(
record.pairname))
valid_record = False
else:
bucket = 'dem-{}s-{}-{}'.format(
args.mode, record.res_str, bp_region.split('-')[0])
custom_path = '/'.join([
path_prefix,
bucket,
res_dir, # e.g. 2m, 50cm, 2m_dsp
record.pairname[:4], # sensor
record.pairname[5:9], # year
record.pairname[9:11], # month
groupid+'.tar' # mode-specific group ID
])
elif args.custom_paths in ('PGC', 'BW'):
# /mnt/pgc/data/elev/dem/setsm/ArcticDEM/region/arcticdem_01_iceland/scenes/
# 2m/WV01_20200630_10200100991E2C00_102001009A862700_2m_v040204/
# WV01_20200630_10200100991E2C00_102001009A862700_
# 504471479080_01_P001_504471481090_01_P001_2_meta.txt
if not region:
logger.error("Pairname not found in region lookup {}, cannot build custom path".format(
record.pairname))
valid_record = False
else:
pretty_project = utils.PROJECTS[region.split('_')[0]]
custom_path = '/'.join([
path_prefix,
pretty_project, # project (e.g. ArcticDEM)
'region',
region, # region
'scenes',
res_dir, # e.g. 2m, 50cm, 2m_dsp
groupid, # strip ID
record.srcfn # file name (meta.txt)
])
elif args.custom_paths == 'CSS':
# /css/nga-dems/setsm/scene/2m/2021/04/21/
# W2W2_20161025_103001005E00BD00_103001005E89F900_2m_v040306
custom_path = '/'.join([
path_prefix,
args.mode, # mode (scene, strip, tile)
res_dir, # e.g. 2m, 50cm, 2m_dsp
record.pairname[:4], # sensor
record.pairname[5:9], # year
record.pairname[9:11], # month
groupid, # mode-specific group ID
record.srcfn # file name (meta.txt)
])
else:
logger.error("Mode {} does not support the specified custom path option,\
skipping record".format(args.mode))
valid_record = False
## Fields for strip DEM
if args.mode == 'strip':
logger.debug(f"Processing strip: {record.stripid}")
attrib_map = {
'DEM_ID': record.stripid,
'STRIPDEMID': record.stripdemid,
'PAIRNAME': record.pairname,
'SENSOR1': record.sensor1,
'SENSOR2': record.sensor2,
'ACQDATE1': record.acqdate1.strftime('%Y-%m-%d'),
'ACQDATE2': record.acqdate2.strftime('%Y-%m-%d'),
'AVGACQTM1': record.avg_acqtime1.strftime("%Y-%m-%d %H:%M:%S") if record.avg_acqtime1 is not None else None,
'AVGACQTM2': record.avg_acqtime2.strftime("%Y-%m-%d %H:%M:%S") if record.avg_acqtime2 is not None else None,
'CATALOGID1': record.catid1,
'CATALOGID2': record.catid2,
'IS_LSF': record.is_lsf,
'IS_XTRACK': record.is_xtrack,
'EDGEMASK': record.mask_tuple[0],
'WATERMASK': record.mask_tuple[1],
'CLOUDMASK': record.mask_tuple[2],
'ALGM_VER': record.algm_version,
'S2S_VER': record.s2s_version,
'RMSE': record.rmse,
'FILESZ_DEM': record.filesz_dem,
'FILESZ_MT': record.filesz_mt,
'FILESZ_OR': record.filesz_or,
'FILESZ_OR2': record.filesz_or2,
'PROJ4': record.proj4,
'EPSG': record.epsg,
'GEOCELL': record.geocell,
}
## Set region
try:
pair_items = pairs[record.pairname]
except KeyError as e:
region = None
bp_region = None
else:
if isinstance(pair_items, str):
region = pair_items
elif isinstance(pair_items, tuple):
region, bp_region = pair_items
else:
logger.error("Pairname-region lookup value cannot be parsed for pairname {}: {}".format(
record.pairname, pair_items))
attrib_map['REGION'] = region
if record.release_version and 'REL_VER' in fld_list:
attrib_map['REL_VER'] = record.release_version
for f, a in utils.field_attrib_map.items():
val = getattr(record, a)
attrib_map[f] = round(val, 6) if val is not None else -9999
## If registration info exists - DEPRECATED
# if args.include_registration:
# if len(record.reginfo_list) > 0:
# for reginfo in record.reginfo_list:
# if reginfo.name == 'ICESat':
# attrib_map["DX"] = reginfo.dx
# attrib_map["DY"] = reginfo.dy
# attrib_map["DZ"] = reginfo.dz
# attrib_map["REG_SRC"] = 'ICESat'
# attrib_map["NUM_GCPS"] = reginfo.num_gcps
# attrib_map["MEANRESZ"] = reginfo.mean_resid_z
## Set path folders for use if path_prefix specified
if path_prefix:
res_dir = record.res_str + '_dsp' if record.is_dsp else record.res_str
if args.custom_paths == 'BP':
# https://blackpearl-data2.pgc.umn.edu/dem-strips-arceua/2m/WV02/2015/05/
# WV02_20150506_1030010041510B00_1030010043050B00_50cm_v040002.tar
if not region:
logger.error("Pairname not found in region lookup {}, cannot build custom path".format(
record.pairname))
valid_record = False
else:
# FIXME: Will we need separate buckets for different s2s version strips (i.e. v4 vs. v4.1)?
bucket = 'dem-{}s-{}'.format(
args.mode, bp_region.split('-')[0])
custom_path = '/'.join([
path_prefix,
bucket,
res_dir, # e.g. 2m, 50cm, 2m_dsp
record.pairname[:4], # sensor
record.pairname[5:9], # year
record.pairname[9:11], # month
groupid + '.tar' # mode-specific group ID
])
elif args.custom_paths in ('PGC', 'BW'):
# /mnt/pgc/data/elev/dem/setsm/ArcticDEM/region/arcticdem_01_iceland/strips_v4/
# 2m/WV01_20200630_10200100991E2C00_102001009A862700_2m_v040204/
# WV01_20200630_10200100991E2C00_102001009A862700_seg1_etc
if not region:
logger.error("Pairname not found in region lookup {}, cannot build custom path".format(
record.pairname))
valid_record = False
else:
pretty_project = utils.PROJECTS[region.split('_')[0]]
custom_path = '/'.join([
path_prefix,
pretty_project, # project (e.g. ArcticDEM)
'region',
region, # region
'strips_v{}'.format(record.s2s_version),
res_dir, # e.g. 2m, 50cm, 2m_dsp
groupid, # strip ID
record.srcfn # file name (meta.txt)
])
elif args.custom_paths == 'CSS':
# /css/nga-dems/setsm/strip/strips_v3/2m/2021/04/21/
# W2W2_20161025_103001005E00BD00_103001005E89F900_2m_v040306
custom_path = '/'.join([
path_prefix,
args.mode, # mode (scene, strip, tile)
'strips_v{}'.format(record.s2s_version),
res_dir, # e.g. 2m, 50cm, 2m_dsp
record.pairname[:4], # sensor
record.pairname[5:9], # year
record.pairname[9:11], # month
groupid, # mode-specific group ID
record.srcfn # file name (meta.txt)
])
else:
logger.error("Mode {} does not support the specified custom path option,\
skipping record".format(args.mode))
valid_record = False
## Fields for tile DEM
if args.mode == 'tile':
logger.debug(f"Processing tile: {record.tileid}")
attrib_map = {
'DEM_ID': record.tileid,
'TILE': record.tile_id_no_res,
'SUPERTILE': record.supertile_id_no_res,
'NUM_COMP': record.num_components,
'FILESZ_DEM': record.filesz_dem,
'EPSG': record.epsg,
}
## Optional attributes
if record.release_version and ('REL_VER' in fld_list or 'RELEASEVER' in fld_list):
attrib_map['REL_VER'] = record.release_version
version = record.release_version
else:
version = 'novers'
attrib_map['DENSITY'] = record.density if record.density is not None else -9999
# if args.include_registration: --DEPRECATED
# if record.reg_src:
# attrib_map["REG_SRC"] = record.reg_src
# attrib_map["NUM_GCPS"] = record.num_gcps
# if record.mean_resid_z:
# attrib_map["MEANRESZ"] = record.mean_resid_z
## Set path folders for use if db_path_prefix specified
if path_prefix:
if args.custom_paths == 'BP':
custom_path = '.'.join([
path_prefix,
record.mode, # mode (scene, strip, tile)
args.project.lower(), # project
record.res, # resolution
version, # version
groupid+'.tar' # mode-specific group ID
])
else:
logger.error("Mode {} does not support the specified custom path option,\
skipping record".format(args.mode))
valid_record = False
## Common fields
if valid_record:
## Common Attributes across all modes
attrib_map['INDEX_DATE'] = datetime.datetime.today().strftime('%Y-%m-%d')
attrib_map['CR_DATE'] = record.creation_date.strftime('%Y-%m-%d')
attrib_map['ND_VALUE'] = record.ndv
if dsp_mode == 'orig':
res = record.dsp_dem_res
else:
res = (record.xres + record.yres) / 2.0
attrib_map['DEM_RES'] = res
## Set location
if path_prefix:
location = custom_path
else:
location = record.srcfp
attrib_map['LOCATION'] = location
## Transform and write geom
src_srs = utils.osr_srs_preserve_axis_order(osr.SpatialReference())
src_srs.ImportFromWkt(record.proj)
if not record.geom:
logger.error('No valid geom found, feature skipped: {}'.format(record.sceneid))
valid_record = False
else:
temp_geom = record.geom.Clone()
transform = osr.CoordinateTransformation(src_srs,tgt_srs)
try:
temp_geom.Transform(transform)
except TypeError as e:
logger.error('Geom transformation failed, feature skipped: {} {}'.format(e, record.sceneid))
valid_record = False
else:
## Get centroid coordinates