<commit_before>#! /usr/bin/env python

#<NAME>
#<EMAIL>

#Utility to co-register two rasters
#This was written in 2012, and needs major cleanup

#Input two uncorrected, unclipped DEMs and a mask for static regions
#Clip both DEMs using mask
#Align the clipped DEMs
#Update georef info and remove elevation offset to DEM2

#Next:
#Need better way to clip input datasets to common mask of static surfaces
#Using gdalwarp is a bit clunky
#This uses PIL to draw polygon
#http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html
#At least remove the temporary clipped files
#mapnik to render shapefile for given extent?
#Define shapefile
#Convert to pbm (or binary ndarray) for specified extent

#To do:
#Clean up path - implement os.path.join(outdir, fn)

import sys
import os
from osgeo import gdal
import numpy as np
import scipy.ndimage

from pygeotools.lib import iolib
from pygeotools.lib import malib
from pygeotools.lib import geolib
from pygeotools.lib import filtlib
from pygeotools.lib import warplib

from demcoreg import coreglib

memdrv = gdal.GetDriverByName('MEM')

def main(argv=None):

    if argv is None:
        argv = sys.argv

    if len(argv) < 3 or len(argv) > 4:
        sys.exit("Usage is %s fixedref_dem source_dem [staticmask.shp]" % argv[0])

    #Input DEMs
    #This is the reference
    dem1_fn = sys.argv[1]
    #This is the source to be moved
    dem2_fn = sys.argv[2]

    print("\nReference: %s" % dem1_fn)
    print("Source: %s" % dem2_fn)

    #Input mask for static areas
    #Note mask must be in same projected coordinate system as input rasters
    #For areas not covering rock, can use the intersection polygon
    mask_fn = None
    if len(sys.argv) == 4:
        mask_fn = sys.argv[3]
        print("Mask: %s\n" % mask_fn)

    dem1_ds = gdal.Open(dem1_fn, gdal.GA_ReadOnly)
    dem2_ds = gdal.Open(dem2_fn, gdal.GA_ReadOnly)

    dem1_gt_orig = np.array(dem1_ds.GetGeoTransform())
    dem2_gt_orig = np.array(dem2_ds.GetGeoTransform())

    #Check to see if input datasets have the same extent (already clipped) and gt
    clip = np.any(dem1_gt_orig != dem2_gt_orig) or np.any(geolib.ds_extent(dem1_ds) != geolib.ds_extent(dem2_ds))

    #outdir = sys.argv[4]
    #mode = sys.argv[3]
    mode = 'ncc'

    dem1_clip_ds = dem1_ds
    dem2_clip_ds = dem2_ds
    
    #These are needed by older filtering code, should be removed (can get from ds)
    dem1_clip_fn = dem1_fn
    dem2_clip_fn = dem2_fn

    if clip:
        if mask_fn is not None:
            #NOTE: this does not resample to ensure both are same res
            #Should be better to clip first, then resample
            print("Clipping input DEMs to input polygon\n")
            dem1_clip_ds = geolib.clip_raster_by_shp(dem1_fn, mask_fn)
            dem2_clip_ds = geolib.clip_raster_by_shp(dem2_fn, mask_fn)
            dem1_clip_fn = dem1_clip_ds.GetFileList()[0]
            dem2_clip_fn = dem2_clip_ds.GetFileList()[0]
            #Need to be careful here, as these are not necessarily set 
            #Also, I think GDAL has trouble when source file for an open ds disappears
            #os.remove(dem1_clip_fn)
            #os.remove(dem2_clip_fn)

        #Make sure the input datasets have the same resolution/extent
        dem1_clip_ds, dem2_clip_ds = warplib.memwarp_multi([dem1_clip_ds, dem2_clip_ds], res='max', extent='intersection')

    dem1_gt = np.array(dem1_clip_ds.GetGeoTransform())
    dem2_gt = np.array(dem2_clip_ds.GetGeoTransform())
    dem1 = iolib.ds_getma(dem1_clip_ds, 1)
    dem2 = iolib.ds_getma(dem2_clip_ds, 1)

    #Compute common mask and apply to all products
    #common_mask = dem1.mask + dem2.mask
    common_mask = malib.common_mask([dem1, dem2]) 

    #Compute difference for unaligned inputs
    diff_euler = np.ma.array(dem1-dem2, mask=common_mask)
    
    print("Elevation difference stats for uncorrected input DEMs")
    diff_stats = malib.print_stats(diff_euler)
    med_bias = diff_stats[5]<commit_msg>fix unused variable in dem_align<commit_after>#! /usr/bin/env python

#<NAME>
#<EMAIL>

#Utility to co-register two rasters
#This was written in 2012, and needs major cleanup

#Input two uncorrected, unclipped DEMs and a mask for static regions
#Clip both DEMs using mask
#Align the clipped DEMs
#Update georef info and remove elevation offset to DEM2

#Next:
#Need better way to clip input datasets to common mask of static surfaces
#Using gdalwarp is a bit clunky
#This uses PIL to draw polygon
#http://geospatialpython.com/2011/02/clip-raster-using-shapefile.html
#At least remove the temporary clipped files
#mapnik to render shapefile for given extent?
#Define shapefile
#Convert to pbm (or binary ndarray) for specified extent

#To do:
#Clean up path - implement os.path.join(outdir, fn)

import sys
import os
from osgeo import gdal
import numpy as np
import scipy.ndimage

from pygeotools.lib import iolib
from pygeotools.lib import malib
from pygeotools.lib import geolib
from pygeotools.lib import filtlib
from pygeotools.lib import warplib

from demcoreg import coreglib

memdrv = gdal.GetDriverByName('MEM')

def main(argv=None):

    if argv is None:
        argv = sys.argv

    if len(argv) < 3 or len(argv) > 4:
        sys.exit("Usage is %s fixedref_dem source_dem [staticmask.shp]" % argv[0])

    #Input DEMs
    #This is the reference
    dem1_fn = sys.argv[1]
    #This is the source to be moved
    dem2_fn = sys.argv[2]

    print("\nReference: %s" % dem1_fn)
    print("Source: %s" % dem2_fn)

    #Input mask for static areas
    #Note mask must be in same projected coordinate system as input rasters
    #For areas not covering rock, can use the intersection polygon
    mask_fn = None
    if len(sys.argv) == 4:
        mask_fn = sys.argv[3]
        print("Mask: %s\n" % mask_fn)

    dem1_ds = gdal.Open(dem1_fn, gdal.GA_ReadOnly)
    dem2_ds = gdal.Open(dem2_fn, gdal.GA_ReadOnly)

    dem1_gt_orig = np.array(dem1_ds.GetGeoTransform())
    dem2_gt_orig = np.array(dem2_ds.GetGeoTransform())

    #Check to see if input datasets have the same extent (already clipped) and gt
    clip = np.any(dem1_gt_orig != dem2_gt_orig) or np.any(geolib.ds_extent(dem1_ds) != geolib.ds_extent(dem2_ds))

    #outdir = sys.argv[4]
    #mode = sys.argv[3]
    mode = 'ncc'

    dem1_clip_ds = dem1_ds
    dem2_clip_ds = dem2_ds
    
    #These are needed by older filtering code, should be removed (can get from ds)
    dem1_clip_fn = dem1_fn
    dem2_clip_fn = dem2_fn

    if clip:
        if mask_fn is not None:
            #NOTE: this does not resample to ensure both are same res
            #Should be better to clip first, then resample
            print("Clipping input DEMs to input polygon\n")
            dem1_clip_ds = geolib.clip_raster_by_shp(dem1_fn, mask_fn)
            dem2_clip_ds = geolib.clip_raster_by_shp(dem2_fn, mask_fn)
            dem1_clip_fn = dem1_clip_ds.GetFileList()[0]
            dem2_clip_fn = dem2_clip_ds.GetFileList()[0]
            #Need to be careful here, as these are not necessarily set 
            #Also, I think GDAL has trouble when source file for an open ds disappears
            #os.remove(dem1_clip_fn)
            #os.remove(dem2_clip_fn)

        #Make sure the input datasets have the same resolution/extent
        dem1_clip_ds, dem2_clip_ds = warplib.memwarp_multi([dem1_clip_ds, dem2_clip_ds], res='max', extent='intersection')

    # dem1_gt = np.array(dem1_clip_ds.GetGeoTransform())
    dem2_gt = np.array(dem2_clip_ds.GetGeoTransform())
    dem1 = iolib.ds_getma(dem1_clip_ds, 1)
    dem2 = iolib.ds_getma(dem2_clip_ds, 1)

    #Compute common mask and apply to all products
    #common_mask = dem1.mask + dem2.mask
    common_mask = malib.common_mask([dem1, dem2]) 

    #Compute difference for unaligned inputs
    diff_euler = np.ma.array(dem1-dem2, mask=common_mask)
    
    print("Elevation difference stats for uncorrected input DEMs")
    diff_stats = malib.print_stats(diff_euler)
    med_bias = diff_stats[5]
{'fraction_non_alphanumeric': 0.08294641736729828, 'fraction_numerical': 0.01932488392521019, 'mean_word_length': 3.345692475463468, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '1032012', 'n_tokens_mistral': 2851, 'n_tokens_neox': 2754, 'n_words': 831}