cffi_geos_transform.py

import numpy as np
from cffi import FFI 
from collections import namedtuple
ffi = FFI()
ffi.cdef("""
enum SCAN_GEOMETRIES {
  GOES,
  GEOS
};

void fgf_to_sat(const double fgf_x, const double fgf_y, const double scale_x, 
                const double offset_x, const double scale_y, const double offset_y, 
                double *lamda, double *theta);
                
void goes_to_geos(double *lamda, double *theta);

void sat_to_earth(const double lamda, const double theta, const double sub_lon_degrees,
                  double *lon_degrees, double *lat_degrees);
                  
void fgf_to_earth_(const double *fgf_x, const double *fgf_y, const double *scale_x, 
                   const double *offset_x, const double *scale_y, const double *offset_y, 
                   const double *sub_lon_degrees, double *lon_degrees, double *lat_degrees);                                  

void earth_to_sat(const double lon_degrees, const double lat_degrees, const double sub_lon_degrees, double *lamda, double *theta);

void sat_to_fgf(const double lamda, const double theta, const double scale_x, const double offset_x, const double scale_y, 
                const double offset_y, double *fgf_x, double *fgf_y);

""")

_so = ffi.dlopen('./libgeos_transform.so')

place = namedtuple('place', ('lon','lat'))

def fgf_to_earth(fgf_x, fgf_y, scale_x, offset_x, scale_y, offset_y, sub_lon_degrees):
    """
    return longitude, latitude    
    """
    lam = ffi.new('double *')
    tht = ffi.new('double *')
    lon = ffi.new('double *')
    lat = ffi.new('double *')
    _so.fgf_to_sat(fgf_x, fgf_y, scale_x, offset_x, scale_y, offset_y, lam, tht)
    print "lam %f tht %f" % (lam[0], tht[0])
    _so.sat_to_earth(lam[0], tht[0], sub_lon_degrees, lon, lat)
    return place(lon[0], lat[0])

DEG2RAD = np.pi / 180.0

def y_x_to_earth(line, column, CFAC, COFF, LFAC, LOFF, sub_lon_degrees):
    'imitate geocoord2pixcoord approach to *FAC/*OFF so we can use info directly from HSD files'
    lambda_ = 2.0**16 * (float(column)-COFF) / CFAC * DEG2RAD
    theta = 2.0**16 * (float(line)-LOFF) / LFAC * DEG2RAD
    lon = ffi.new('double *')
    lat = ffi.new('double *')
    _so.sat_to_earth(lambda_, theta, sub_lon_degrees, lon, lat)
    return place(lon[0], lat[0])


def hsd_test(line=1856, column=1858):
    sub_lon = 140.7
    CFAC = 20466275
    LFAC = 20466275
    COFF = 2750.5
    LOFF = 2750.5
    return lc_to_earth(line, column, CFAC, COFF, LFAC, LOFF, sub_lon)




class fgf(object):
    _xoyo = None
    
    def __init__(self, cfac, lfac, coff, loff, sublon):
        self._xoyo = xoyo = ffi.new('double[5]')
        xoyo[0] = cfac
        xoyo[1] = coff
        xoyo[2] = lfac
        xoyo[3] = loff
        xoyo[4] = sublon
    
    def __call__(self, fgf_x, fgf_y):
        lam = ffi.new('double *')
        tht = ffi.new('double *')
        lon = ffi.new('double *')
        lat = ffi.new('double *')
        _so.fgf_to_sat(fgf_x, fgf_y, self._xoyo[0], self._xoyo[1], self._xoyo[2], self._xoyo[3], lam, tht)
        _so.sat_to_earth(lam[0], tht[0], self._xoyo[4], lon, lat)
        return lon[0], lat[0]