Commit 3daad0eb authored by R.K.Garcia's avatar R.K.Garcia
Browse files

fix package grammar after getting my knuckles rapped by the teacher T_T

parent 9e52c7be
HIMAWARI_INTERFACE = """# 1 "../include/HimawariTypes.h"
# 1 "<built-in>" 1
# 1 "<built-in>" 3
# 325 "<built-in>" 3
# 1 "<command line>" 1
# 1 "<built-in>" 2
# 1 "../include/HimawariTypes.h" 2
typedef void (*HimawariLogFunction)(int severity, const char *explanation);
struct HimawariSentinels {
uint16_t count_invalid;
uint16_t count_outside_scan;
double derived_invalid;
double derived_outside_scan;
};
struct HimawariCoordsScale {
double my;
double by;
double mx;
double bx;
};
struct HimawariTime {
int32_t year,
month,
day,
hour,
minute,
second,
microsecond;
int32_t doy;
double unix_time;
double mjd;
};
enum HimawariBandType {
AHI_BAND_UNKNOWN = 0,
AHI_BAND_ALBEDO = 1,
AHI_BAND_BRIGHTNESS_TEMP = 2
};
struct HimawariNavigation {
double sub_lon;
int32_t LFAC, CFAC;
float LOFF, COFF;
double distance_from_earth_center_to_virtual_satellite;
double earth_equatorial_radius;
double earth_polar_radius;
float ssp_latitude, ssp_longitude;
double nadir_latitude, nadir_longitude;
double distance_from_earth_center_to_satellite;
};
struct HimawariCalibration {
double rad_m;
double rad_b;
double bt_c0_or_albedo_cprime;
double bt_c1;
double bt_c2;
double rad_c0;
double rad_c1;
double rad_c2;
double c;
double h;
double k;
};
struct HimawariMetadata {
int32_t lines, columns;
int32_t lines_res_meters,columns_res_meters;
int32_t observation_timeline;
struct HimawariTime start_time, end_time;
int32_t begin_line, end_line;
int32_t available_lines;
int32_t band;
int32_t band_type;
double central_wavelength;
int resolution;
int32_t wmo_sat_id;
char satellite_name[17];
};
struct HimawariSourceLocation {
int32_t line_offset, column_offset;
int32_t line_stride, column_stride;
int32_t lines, columns;
int32_t nav_options;
};
struct HimawariDestinationLocation {
int32_t line_increment, column_increment;
int32_t downsample_to_match_band;
};
extern const struct HimawariSentinels DEFAULT_HIMAWARI_SENTINELS;
extern const struct HimawariSentinels DEFAULT_HCAST_SENTINELS;
# 1 "../include/HimawariConstants.h"
# 1 "<built-in>" 1
# 1 "<built-in>" 3
# 325 "<built-in>" 3
# 1 "<command line>" 1
# 1 "<built-in>" 2
# 1 "../include/HimawariConstants.h" 2
const int LOG_WARNING=0,
LOG_ERROR=-1,
LOG_CRITICAL=-2,
LOG_INFO=1,
LOG_DEBUG=2;
const int HIMAWARI_RES_0P5KM = 500,
HIMAWARI_RES_1P0KM = 1000,
HIMAWARI_RES_2P0KM = 2000,
HIMAWARI_RES_4P0KM = 4000
;
const int NAVOPT_CORRECT_DATA = 1;
const int NAVOPT_CORRECT_FGF = 2;
const int NAVOPT_DEFAULT = NAVOPT_CORRECT_DATA;
# 1 "../include/HimawariErrors.h"
# 1 "<built-in>" 1
# 1 "<built-in>" 3
# 325 "<built-in>" 3
# 1 "<command line>" 1
# 1 "<built-in>" 2
# 1 "../include/HimawariErrors.h" 2
const int ERROR_UNIMPLEMENTED = -9999;
const int ERROR_NOT_AVAILABLE = -999;
const int ERROR_ATTACH_DATA_BLOCK=-2,
ERROR_ATTACH_PROJ_BLOCK=-3,
ERROR_ATTACH_NAV_BLOCK=-4,
ERROR_ATTACH_CAL_BLOCK=-5,
ERROR_ATTACH_INTERCAL_BLOCK=-6,
ERROR_ATTACH_SEGMENT_BLOCK=-7,
ERROR_ATTACH_NAV_CORR_BLOCK=-8,
ERROR_ATTACH_OBS_TIME_BLOCK=-9,
ERROR_ATTACH_ERROR_BLOCK=-10;
const int ERROR_INVALID_SCENE = -100000;
const int ERROR_TIME_CONVERSION = -100001;
const int ERROR_UNSUPPORTED_STRIDE = -100002;
const int ERROR_COPY_NOT_ALBEDO=-11,
ERROR_COPY_NOT_BRIGHTNESS_TEMP=-12,
ERROR_COPY_OUTSIDE_BOUNDS=-13;
const int WARNING_SOME_DATA_MISSING = 1;
# 1 "../include/HimawariFunctions.h"
# 1 "<built-in>" 1
# 1 "<built-in>" 3
# 325 "<built-in>" 3
# 1 "<command line>" 1
# 1 "<built-in>" 2
# 1 "../include/HimawariFunctions.h" 2
struct HimawariScene *openHimawariScene(
const char *path_stem,
HimawariLogFunction log_function,
const struct HimawariSentinels *sentinels
);
# 24 "../include/HimawariFunctions.h"
int queryHimawariSceneMetadata(const struct HimawariScene *scene,
struct HimawariMetadata *metadata)
;
int queryHimawariSceneNavigation(const struct HimawariScene *scene,
struct HimawariNavigation *navigation
);
int queryHimawariSceneCalibration(const struct HimawariScene *scene,
struct HimawariCalibration *cal
);
int queryHimawariSceneLineTimeOffsets(const struct HimawariScene *scene,
double *seconds
);
# 68 "../include/HimawariFunctions.h"
int copyHimawariSceneCounts(const struct HimawariScene *it,
const struct HimawariSourceLocation *from_where,
const struct HimawariDestinationLocation *to_where,
float *output
);
int copyHimawariSceneRadiances(const struct HimawariScene *it,
const struct HimawariSourceLocation *from_where,
const struct HimawariDestinationLocation *to_where,
float *output
);
int copyHimawariSceneBrightnessTemps(const struct HimawariScene *it,
const struct HimawariSourceLocation *from_where,
const struct HimawariDestinationLocation *to_where,
float *output
);
int copyHimawariSceneAlbedos(const struct HimawariScene *it,
const struct HimawariSourceLocation *from_where,
const struct HimawariDestinationLocation *to_where,
float *output
);
# 115 "../include/HimawariFunctions.h"
int copyHimawariSceneCoords(const struct HimawariScene *it,
const struct HimawariSourceLocation *from_where,
const struct HimawariDestinationLocation *to_where,
float *y_output,
float *x_output,
struct HimawariCoordsScale *scale
);
int closeHimawariScene(struct HimawariScene *it);
int averageHimawariFloatArray(
int32_t factor,
float sentinel,
const float *source,
int32_t source_lines, int32_t source_columns,
int32_t source_line_increment, int32_t source_column_increment,
int32_t source_first_line, int32_t source_first_column,
int32_t source_band,
float *dest,
int32_t dest_lines, int32_t dest_columns,
int32_t dest_line_increment, int32_t dest_column_increment,
int32_t dest_band
);
"""
This diff is collapsed.
This diff is collapsed.
This diff is collapsed.
#!/usr/bin/env python
"""
Check whether a given set of bands and sections is available, returning suitable success/failure result code
# satisfied requirements
python gate.py 'HS_H08_20130710_0300_B{bb}_FLDK_R??_S{ss}10.DAT' '[1,2,3,4,5]' '[1,2,3]' && echo hooray
hooray
# unsatisfied requirements
python gate.py 'HS_H08_20130710_0300_B{bb}_FLDK_R??_S{ss}10.DAT' '[1,2,3,4,5,9]' '[1,2,3,11]' || echo nuuu
nuuu
"""
import sys, os
from glob import glob
def gate(stem, bands, sections):
"""
check that all desired bands and segments are available, return True if so
>>> gate( '/path/to/scene/HS_H08_yyyymmdd_hhmm_B{bb}_FLDK_R??_S{ss}10.dat', [1,2,3,4,5], [1,2,3])
True
"""
bands = ['%02d' % x for x in bands]
print(bands)
sections = ['%02d' % x for x in sections]
print (sections)
exists = lambda pat: len(list(glob(pat)))>0
for permutation in [{'bb':bb, 'ss':ss} for bb in bands for ss in sections]:
path = stem.format(**permutation)
if not exists(path):
print "missing {}".format(path)
return False
else:
print "found {}".format(path)
return True
if __name__=='__main__':
assert(len(sys.argv)==4)
stem = sys.argv[1]
bands = eval(sys.argv[2])
sections = eval(sys.argv[3])
sys.exit(0 if gate(stem, bands, sections) else 1)
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]
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
cgms_nav
~~~~~~~~
Module to navigate image coordinates to Fixed-Grid lat/lon earth coordinates
reference: CGMS 03 (Coordination Group for Meteorological Satellites) LRIT/HRIT Global Specification
reference: pixcoord2geocoord.pro IDL subroutine, wstraka@ssec.wisc.edu
reference: geos_transform.c, wstraka@ssec.wisc.edu / graemem@ssec.wisc.edu
:copyright: 2014 by University of Wisconsin Regents, see AUTHORS for more details
:license: GPLv3, see LICENSE for more details
"""
__author__ = 'rayg'
__docformat__ = 'reStructuredText'
import os, sys
import logging, unittest, argparse
import numpy as np
from collections import namedtuple
LOG = logging.getLogger(__name__)
geocoords = namedtuple('geocoord', ('latitude', 'longitude', 'y', 'x'))
DEG2RAD = np.pi / 180.0
def hsd_lcof2yx(lc, off, fac):
"convert HSD grid line/column to y/x"
return (float(2**16) * (float(lc) - off)) / float(fac) * DEG2RAD
# no conversion, line is fgf_y and column is fgf_x
yx2yx = lambda q,off,fac: q
def cgms_nav(line, column, coff, loff, cfac, lfac, sub_lon,
sat_height=42164.0, r_eq=6378.1370, r_pol=6356.7523, lcof2yx=None):
"""
perform line/column to lat/lon for array inputs
:param line: array of line indices, can be fractional
:param column: array of column indices, can be fractional
:param coff: column offset constant
:param loff: line offset constant
:param cfac: column factor constant
:param lfac: line factor constant
:param sub_lon: satellite effective subpoint for this projection, degrees
:param sat_height: satellite height relative to center of geoid, km
:param r_eq: equatorial radius of earth geoid, kilometers
:param r_pol: polar radius of earth geoid, kilometers
:param lcof2yx: function converting (line or column), offset, factor => FGF y or x angle in radians
:return: geocoords tuple including masked arrays of latitude and longitude
adapted from WCS3 pixcoord2geocoord IDL function.
"""
LOG.debug(locals())
# calculate viewing angle of the satellite by use of the equation
# on page 28, Ref [1].
if lcof2yx is None:
y = line
x = column
else:
y = lcof2yx(line, loff, lfac)
x = lcof2yx(column, coff, cfac)
del line, column
# Only added to speed things up---- WCS3
cos_x_cos_y = np.cos(x) * np.cos(y)
cos_y_cos_y = np.power(np.cos(y), 2)
sin_y_sin_y = np.power(np.sin(y), 2)
sin_x_cos_y = np.sin(x) * np.cos(y)
# now calculate the inverse projection using equations on page 25, Ref. [1]
# first check for visibility, whether the pixel is located on the Earth
# surface or in space.
# To do this calculate the argument to sqrt of "sd", which is named "sa".
# If it is negative then the sqrt will return NaN and the pixel will be
# located in space, otherwise all is fine and the pixel is located on the
# Earth surface.
sa = np.power(cos_x_cos_y * sat_height, 2) - (cos_y_cos_y + 1.006803 * sin_y_sin_y) * 1737121856.0
# take care if the pixel is in space, that an error code will be returned
mask = sa < 0.0
del sa
cos_x_cos_y = np.ma.masked_array(cos_x_cos_y, mask, fill_value=np.nan)
cos_y_cos_y = np.ma.masked_array(cos_y_cos_y, mask, fill_value=np.nan)
sin_y_sin_y = np.ma.masked_array(sin_y_sin_y, mask, fill_value=np.nan)
sin_x_cos_y = np.ma.masked_array(sin_x_cos_y, mask, fill_value=np.nan)
sin_y = np.ma.masked_array(np.sin(y), mask, fill_value=np.nan)
del mask
# now calculate the rest of the formulas using eq. on page 25 Ref [1]
sd = np.sqrt(np.power(cos_x_cos_y * sat_height, 2) - (cos_y_cos_y + 1.006803 * sin_y_sin_y) * 1737121856.)
sn = (cos_x_cos_y * sat_height - sd) / (cos_y_cos_y + 1.006803 * sin_y_sin_y)
del sd
s1 = -sn * cos_x_cos_y + sat_height
s2 = sn * sin_x_cos_y
s3 = -sn * sin_y
del sin_y
sxy = np.sqrt( np.power(s1, 2) + np.power(s2, 2) )
# using the previous calculations now the inverse projection can be
# calculated, which means calculating the lat./long. from the pixel
# row and column by equations on page 25, Ref [1].
longi = np.arctan2(s2, s1)
lati = np.arctan2((1.006803*s3), sxy)
del s1, s2, s3, sxy
# convert from radians into degrees
latitude = lati/np.pi*180.0
longitude = (longi/np.pi*180.0) + sub_lon
return geocoords(latitude=latitude, longitude=longitude, y=y, x=x)
def cgms_proj4(sub_lon_deg, sat_height_km, r_eq_km, r_pol_km, x_0 = 0, y_0 = 0):
"""
Convert CGMS03 to Proj.4 projection string
:param sub_lon_deg: effective sub-satellite longitude for the projected data
:param sat_height_km: effective satellite height in km
:param r_eq_km: equatorial geoid radius in km
:param r_pol_km: polar geoid radius in km
:param x_0: false easting in km
:param y_0: false northing in km
:return: proj4 string
"""
sat_height = sat_height_km * 1000.0
a = r_eq_km * 1000.0
b = r_pol_km * 1000.0
x_0 *= 1000.0
y_0 *= 1000.0
return "+proj=geos +lon_0={sub_lon_deg:.4f} +h={sat_height:.5f} +x_0={x_0:d} +y_0={y_0:d} +a={a:.4f} +b={b:.4f} +units=m".format(globals())
import os,sys
import HimawariScene as hsd
D = {}
for fn in sys.argv[1:]:
hs = hsd.HimawariScene(fn)
meta = hs.metadata
nav = hs.navigation
q = dict(CFAC=nav.CFAC, LFAC=nav.LFAC, COFF=nav.COFF, LOFF=nav.LOFF)
# D[(meta.band, meta.lines)] = q