1# This file is part of the Astrometry.net suite. 2# Licensed under a 3-clause BSD style license - see LICENSE 3from __future__ import print_function 4import matplotlib 5matplotlib.use('Agg') 6 7import sys 8from optparse import * 9 10import numpy as np 11from pylab import * 12from numpy import * 13#from astrometry.util.sip import * 14from astrometry.util.util import * 15 16def plot_distortions(wcsfn, ex=1, ngridx=10, ngridy=10, stepx=10, stepy=10): 17 wcs = Sip(wcsfn) 18 W,H = wcs.wcstan.imagew, wcs.wcstan.imageh 19 20 xgrid = np.linspace(0, W, ngridx) 21 ygrid = np.linspace(0, H, ngridy) 22 X = np.linspace(0, W, int(ceil(W/stepx))) 23 Y = np.linspace(0, H, int(ceil(H/stepy))) 24 25 xlo,xhi,ylo,yhi = 0,W,0,H 26 27 for x in xgrid: 28 DX,DY = [],[] 29 xx,yy = [],[] 30 UX,UY = [],[] 31 for y in Y: 32 dx,dy = wcs.get_distortion(x, y) 33 ux,uy = wcs.get_undistortion(dx, dy) 34 print('x,y', (x,y), 'dx,dy', (dx,dy), 'ux,uy', (ux,uy)) 35 xx.append(x) 36 yy.append(y) 37 DX.append(dx) 38 DY.append(dy) 39 UX.append(ux) 40 UY.append(uy) 41 DX = np.array(DX) 42 DY = np.array(DY) 43 UX = np.array(UX) 44 UY = np.array(UY) 45 xx = np.array(xx) 46 yy = np.array(yy) 47 EX = DX + ex * (DX - xx) 48 EY = DY + ex * (DY - yy) 49 plot(xx, yy, 'k-', alpha=0.5) 50 plot(EX, EY, 'r-') 51 plot(UX, UY, 'b-', alpha=0.5) 52 xlo = min(xlo, min(EX)) 53 xhi = max(xhi, max(EX)) 54 ylo = min(ylo, min(EY)) 55 yhi = max(yhi, max(EY)) 56 57 for y in ygrid: 58 DX,DY = [],[] 59 xx,yy = [],[] 60 UX,UY = [],[] 61 for x in X: 62 dx,dy = wcs.get_distortion(x, y) 63 ux,uy = wcs.get_undistortion(dx, dy) 64 DX.append(dx) 65 DY.append(dy) 66 UX.append(ux) 67 UY.append(uy) 68 xx.append(x) 69 yy.append(y) 70 DX = np.array(DX) 71 DY = np.array(DY) 72 xx = np.array(xx) 73 yy = np.array(yy) 74 EX = DX + ex * (DX - xx) 75 EY = DY + ex * (DY - yy) 76 plot(xx, yy, 'k-', alpha=0.5) 77 plot(EX, EY, 'r-') 78 plot(UX, UY, 'b-', alpha=0.5) 79 xlo = min(xlo, min(EX)) 80 xhi = max(xhi, max(EX)) 81 ylo = min(ylo, min(EY)) 82 yhi = max(yhi, max(EY)) 83 84 plot([wcs.wcstan.crpix[0]], [wcs.wcstan.crpix[1]], 'rx') 85 86 #axis([0, W, 0, H]) 87 axis('scaled') 88 axis([xlo,xhi,ylo,yhi]) 89 #axis('tight') 90 91if __name__ == '__main__': 92 parser = OptionParser(usage='%prog [options] <wcs-filename> <plot-filename>') 93 parser.add_option('-e', '--ex', '--exaggerate', dest='ex', type='float', help='Exaggerate the distortion by this factor') 94 #parser.add_option('-s', '--scale', dest='scale', type='float', help='Scale the 95 parser.add_option('-n', dest='nsteps', type='int', help='Number of grid lines to plot') 96 97 parser.set_defaults(ex=1.) 98 99 opt,args = parser.parse_args() 100 101 if len(args) != 2: 102 parser.print_help() 103 sys.exit(-1) 104 105 wcsfn = args[0] 106 outfn = args[1] 107 108 args = {} 109 if opt.ex is not None: 110 args['ex'] = opt.ex 111 if opt.nsteps is not None: 112 args['ngridx'] = opt.nsteps 113 args['ngridy'] = opt.nsteps 114 115 clf() 116 plot_distortions(wcsfn, **args) 117 tt = 'SIP distortions: %s' % wcsfn 118 if opt.ex != 1: 119 tt += ' (exaggerated by %g)' % opt.ex 120 title(tt) 121 savefig(outfn) 122 123