ISCE_RTC_geocode_Test.py

Helen Baldwin, 05/04/2020 08:48 pm

Download (3.9 kB)

 
1
#!/usr/bin/env python2
2
# -*- coding: utf-8 -*-
3
"""
4
Created on Thu Apr 30 14:21:16 2020
5

6
@author: fsh
7
"""
8

    
9
from osgeo import gdal, ogr, osr
10
import matplotlib.pyplot as plt
11

    
12
gdal.UseExceptions()
13

    
14
#Utility function to load data
15
def loadData(infile, band=1):
16
    ds = gdal.Open(infile, gdal.GA_ReadOnly)
17
    #Data array
18
    data = ds.GetRasterBand(band).ReadAsArray
19
    #Map extent
20
    trans = ds.GetGeoTransform()
21
    xsize = ds.RasterXSize
22
    ysize = ds.RasterYSize
23
    extent = [trans[0], trans[0], + xsize * trans[1], trans[3], 
24
              ysize * trans[5], trans[3]]
25
    
26
    ds = None
27
    return data, extent
28

    
29

    
30
#gdalwarp [--help-general] [--formats]
31
#    [-s_srs srs_def] [-t_srs srs_def] [-to "NAME=VALUE"]* [-novshiftgrid]
32
#    [-order n | -tps | -rpc | -geoloc] [-et err_threshold]
33
#    [-refine_gcps tolerance [minimum_gcps]]
34
#    [-te xmin ymin xmax ymax] [-te_srs srs_def]
35
#    [-tr xres yres] [-tap] [-ts width height]
36
#    [-ovr level|AUTO|AUTO-n|NONE] [-wo "NAME=VALUE"] [-ot Byte/Int16/...] [-wt Byte/Int16]
37
#    [-srcnodata "value [value...]"] [-dstnodata "value [value...]"]
38
#    [-srcalpha|-nosrcalpha] [-dstalpha]
39
#    [-r resampling_method] [-wm memory_in_mb] [-multi] [-q]
40
#    [-cutline datasource] [-cl layer] [-cwhere expression]
41
#    [-csql statement] [-cblend dist_in_pixels] [-crop_to_cutline]
42
#    [-of format] [-co "NAME=VALUE"]* [-overwrite]
43
#    [-nomd] [-cvmd meta_conflict_value] [-setci] [-oo NAME=VALUE]*
44
#    [-doo NAME=VALUE]*
45
#    srcfile* dstfile
46

    
47
def geocodeUsingGdalWarp(infile, latfile, lonfile, outfile,
48
                         insrs=4326, outsrs=None,
49
                         spacing=None, fmt='GTiff', bounds=None,
50
                         method='near'):
51
    '''
52
    Geocode a swath file using corresponding lat, lon files
53
    '''
54
    sourcexmltmpl = '''    <SimpleSource>
55
      <SourceFilename>{0}</SourceFilename>
56
      <SourceBand>{1}</SourceBand>
57
    </SimpleSource>'''
58
    
59
    driver = gdal.GetDriverByName('VRT')
60
    tempvrtname = '/media/sf_share/RTC_4-15-2020/mdd/Geocode_test1/geocode.vrt'
61
    inds = gdal.OpenShared(infile, gdal.GA_ReadOnly)
62
    
63
    tempds = driver.Create(tempvrtname, inds.RasterXSize, inds.RasterYSize, 0)
64
    
65
    for ii in range(inds.RasterCount):
66
        band = inds.GetRasterBand(1)
67
        tempds.AddBand(band.DataType)
68
        tempds.GetRasterBand(ii+1).SetMetadata({'source_0': sourcexmltmpl.format(infile, ii+1)}, 'vrt_sources')
69
  
70
    sref = osr.SpatialReference()
71
    sref.ImportFromEPSG(insrs)
72
    srswkt = sref.ExportToWkt()
73

    
74
    tempds.SetMetadata({'SRS' : srswkt,
75
                        'X_DATASET': lonfile,
76
                        'X_BAND' : '1',
77
                        'Y_DATASET': latfile,
78
                        'Y_BAND' : '1',
79
                        'PIXEL_OFFSET' : '0',
80
                        'LINE_OFFSET' : '0',
81
                        'PIXEL_STEP' : '1',
82
                        'LINE_STEP' : '1'}, 'GEOLOCATION')
83
    
84
    band = None
85
    tempds = None 
86
    inds = None
87
    
88
    if spacing is None:
89
        spacing = [None, None]
90
    
91
    warpOptions = gdal.WarpOptions(format=fmt,
92
                                   xRes=spacing[0], yRes=spacing[0],
93
                                   dstSRS=outsrs, outputBounds = bounds, 
94
                                   resampleAlg=method, geoloc=True)
95
    gdal.Warp(outfile, tempvrtname, options=warpOptions)
96

    
97
###Create and visualize geocoded output
98
geocodeUsingGdalWarp('/media/sf_share/RTC_4-15-2020/mdd/output/gamma_VH.10alks_10rlks.img.vrt',
99
                     '/media/sf_share/RTC_4-15-2020/mdd/geometry/lat.rdr.vrt',
100
                     '/media/sf_share/RTC_4-15-2020/mdd/geometry/lon.rdr.vrt',
101
                     '/media/sf_share/RTC_4-15-2020/mdd/Geocode_test1/geocoded_mdd-rtc.tif',
102
                     spacing=[0.0005, 0.0005])
103

    
104
#geodata, geoext = loadData('/media/sf_share/RTC_4-15-2020/mdd/Geocode_test1/geocoded_mdd-rtc.tif', band=2)
105

    
106
#plt.figure()
107
#plt.imshow(geodata, extent=geoext, cmap='gray')
108
#plt.show()
109