These are chat archives for BernieenJos/weetikveel

23rd
Apr 2015
JosvanderSpek
@JosvanderSpek
Apr 23 2015 15:35
Hoi
iambernie
@iambernie
Apr 23 2015 15:35
hio
JosvanderSpek
@JosvanderSpek
Apr 23 2015 15:36

import argparse
import pyfits
import numpy
import sys
import os
import shutil as sh
import subprocess as sub

def get_arguments():
parser = argparse.ArgumentParser()
parser.add_argument("fitsfiles", nargs = "+")
parser.add_argument("--grid", "-g", type = int, nargs = 2, required = True)
parser.add_argument("--outputdir","-o")
args = parser.parse_args()
return args

def get_boundaries(pixels, dim):
xs = [0]
gridsize = pixels // dim
while xs[-1] != pixels:
if (xs[-1] + gridsize) >= pixels :
xs.append(pixels)
else :
xs.append(xs[-1] + gridsize)
return xs

args = get_arguments()
all_segments = []
dimx = args.grid[0]
dimy = args.grid[1]

def make_segmenten(fitsfile, xs, ys):
img = pyfits.open(fitsfile)[0].data
for i in range(len(xs) - 1):
for j in range(len(ys) - 1):
segment = img[xs[i]:xs[i+1], ys[j]:ys[j+1]]
(pyfits.PrimaryHDU(segment)).writeto(abspathdirname+"/segment_%i%i"%(j+1,i+1)+".fits", clobber = True)

for fits in args.fitsfiles:
naam = os.path.basename(fits)
name_dir = os.path.splitext(naam)[0]
abspath_dirname = args.outputdir+"/"+name_dir
if not os.path.exists(abspath_dirname):
os.makedirs(abspath_dirname)
header = pyfits.open(fits)[0].header
pixelsx = header['NAXIS2']
pixelsy = header['NAXIS1']
xs = get_boundaries(pixelsx, dimx)
ys = get_boundaries(pixelsy, dimy)
make_segmenten(fits, xs, ys)

for i in range(len(xs) - 1):
for j in range(len(ys) - 1):
sh.move("segments/jos1/segment%i%i"%(j+1,i+1)+".fits", os.getcwd() + "/infile0.fits")
sh.move("segments/jos2/segment%i%i"%(j+1,i+1)+".fits", os.getcwd() + "/infile1.fits")
os.chdir('../register2')
sub.call('interp.csh')
sub.call(['../bin/mrj_phot', '../images2/interp_infile0.fits', '../images2/interp_infile1.fits']) # segment_i_j.fits is nu verwerkt tot register2/conv.fits
sh.move('conv.fits', "processedsegments/segment_%i%i"%(j+1,i+1)+".fits")
os.chdir('../images2')

Als split.py gerund wordt in images2, dan plaatst hij de segmenten in images2/Segments/jos/segment__*

Voordat interp gerund kan worden, moeten de te runnen segmenten als infile*.fits in images2 staan.

Dit kan wellicht met een nieuw programma dat images2/Segments/jos1/segment_1_1.fits en images2/Segments/jos2/segment_1_1.fits

kopieert/verplaatst naar images2/infile0.fits respectievelijk images2/infile1.fits, vervolgens ../register2/interp.csh runt^

(hetgeen de files images2/interp_infile0.fits en images2/interp_infile1.fits genereert) en vervolgens

../bin/mrj_phot ../images2/interp_infile0.fits ../images2/interp_infile1.fits runt

^interp.csh kan niet vanuit images2 gerund worden! (i.e. in images2 kan men niet zeggen ../register2/interp.csh)

Dus het programma moet de volgende dingen doen:

run split.py in images2. Resultaat: images2/Segments/jos/segment__*

kopieer/verplaats images2/Segments/jos1/segment_1_1.fits en images2/Segments/jos2/segment_1_1.fits ## Uiteraard is dit de eerste iteratie in de uiteindelijke for loop

naar images2/infile0.fits respectievelijk images2/infile1.fits

os.chdir('../register2')

sub.call('interp.csh') Resultaat: images2/interp_infile0.fits en images2/interp_infile1.fits.

sub.call('../bin/mrj_phot') ../images2/interp_infile0.fits ../images2/interp_infile1.fits (ARGUMENTEN ???) ## Op de een of andere manier; sub.call('../bin/mrj_phot')

## werkt namelijk niet, en we kunnen ook niet zomaar

## naar ../bin want dan vindt hij process_config weer niet.

Anyway, stel dat dit gelukt zou zijn, dan zou het resultaat de afbeelding register2/conv.fits creeren.

Kopieer/verplaats deze afbeelding naar een map register2/processed_segments/segment_1_1.fits

De tweede iteratie van de for loop kan segment_1_2.fits gaan behandelen, etc.

Als alle iteraties geweest zijn, en alle verwerkte afbeeldingen in de map re

'''
import argparse
import pyfits
import numpy
import sys
import os
import shutil as sh
import subprocess as sub

def get_arguments():
    parser = argparse.ArgumentParser()
    parser.add_argument("fitsfiles", nargs = "+")
    parser.add_argument("--grid", "-g", type = int, nargs = 2, required = True)
    parser.add_argument("--outputdir","-o")
    args = parser.parse_args()
    return args


def get_boundaries(pixels, dim):
    xs  = [0]
    gridsize = pixels // dim
    while xs[-1] != pixels: 
        if (xs[-1] + gridsize) >= pixels :
            xs.append(pixels)
        else :  
            xs.append(xs[-1] + gridsize)
    return xs

args = get_arguments()
all_segments = []
dimx    = args.grid[0] 
dimy    = args.grid[1]


def make_segmenten(fitsfile, xs, ys):
    img = pyfits.open(fitsfile)[0].data
    for i in range(len(xs) - 1):
        for j in range(len(ys) - 1):
            segment = img[xs[i]:xs[i+1], ys[j]:ys[j+1]] 
            (pyfits.PrimaryHDU(segment)).writeto(abspath_dirname+"/segment_%i_%i"%(j+1,i+1)+".fits", clobber = True)


for fits in args.fitsfiles:
    naam    = os.path.basename(fits)
    name_dir = os.path.splitext(naam)[0]
    abspath_dirname = args.outputdir+"/"+name_dir
    if not os.path.exists(abspath_dirname):
        os.makedirs(abspath_dirname)
    header  = pyfits.open(fits)[0].header
    pixelsx = header['NAXIS2']
    pixelsy = header['NAXIS1']
    xs = get_boundaries(pixelsx, dimx)
    ys = get_boundaries(pixelsy, dimy)
    make_segmenten(fits, xs, ys)

for i in range(len(xs) - 1):
    for j in range(len(ys) - 1):
        sh.move("segments/jos1/segment_%i_%i"%(j+1,i+1)+".fits", os.getcwd() + "/infile0.fits")
        sh.move("segments/jos2/segment_%i_%i"%(j+1,i+1)+".fits", os.getcwd() + "/infile1.fits")
        os.chdir('../register2')
        sub.call('interp.csh')    
        sub.call(['../bin/mrj_phot', '../images2/interp_infile0.fits', '../images2/interp_infile1.fits']) # segment_i_j.fits is nu verwerkt tot register2/conv.fits
        sh.move('conv.fits', "processed_segments/segment_%i_%i"%(j+1,i+1)+".fits")
        os.chdir('../images2')
iambernie
@iambernie
Apr 23 2015 16:09

Voorbeeld hstack/vstack:

>>> a = numpy.arange(16).reshape(4,4)
>>> a
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])
>>> s1 = a[0:2,0:2]
>>> s2 = a[0:2,2:4]
>>> s3 = a[2:4,0:2]
>>> s4 = a[2:4,2:4]
>>> s1
array([[0, 1],
       [4, 5]])
>>> s2
array([[2, 3],
       [6, 7]])
>>> s3
array([[ 8,  9],
       [12, 13]])

>>> numpy.vstack([s1,s2])
array([[0, 1],
       [4, 5],
       [2, 3],
       [6, 7]])

>>> numpy.hstack([s1,s2])
array([[0, 1, 2, 3],
       [4, 5, 6, 7]])

>>> numpy.hstack([s3,s4])
array([[ 8,  9, 10, 11],
       [12, 13, 14, 15]])

>>> numpy.vstack( [numpy.hstack([s1,s2]), numpy.hstack([s3,s4]) ])
array([[ 0,  1,  2,  3],
       [ 4,  5,  6,  7],
       [ 8,  9, 10, 11],
       [12, 13, 14, 15]])

Een mogelijke oplossing zou bijv iets van een (evt. dubbele) for-loop te schrijven die iteratief hstack en vstack aanroept.

iambernie
@iambernie
Apr 23 2015 16:35
En zo kun je een header keyword toevoegen aan je segmenten:
>>> a = numpy.arange(16).reshape(4,4)
>>> s1 = a[0:2,0:2]
>>> s2 = a[0:2,2:4]
>>> s3 = a[2:4,0:2]
>>> s4 = a[2:4,2:4]
>>> pHDU1 = pyfits.PrimaryHDU(s1)
>>> pHDU1.header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                   64 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                    2                                                  
NAXIS2  =                    2                                                  
EXTEND  =                    T                                                  
>>> pHDU1.header['SEGMENTX'] = 1
>>> pHDU1.header['SEGMENTY'] = 1
>>> pHDU1.header
SIMPLE  =                    T / conforms to FITS standard                      
BITPIX  =                   64 / array data type                                
NAXIS   =                    2 / number of array dimensions                     
NAXIS1  =                    2                                                  
NAXIS2  =                    2                                                  
EXTEND  =                    T                                                  
SEGMENTX=                    1                                                  
SEGMENTY=                    1                                                  
>>> pHDU1.writeto('seg1.fits')
>>>
Ik zou voor de zekerheid ook header KEYWORDS beperken tot 8 karakters.