temp commit

ancillary.pyx

+# cython: language_level=3
+# cython: c_string_Type=unicode, c_string_encoding=utf8
+
+cdef extern void get_Reynolds_SST(float *, float *, char *, char *, float *)
+cdef extern void get_NDVI_background(float *, float *, char *, char *, float *)
+cdef extern void get_Olson_eco(float *, float *, char *, unsigned char *)
+cdef extern void get_GEOS(float *, float *, char *, char *, char *, char *, char *, char *, char *,
+                          float *, float *, float *, float *, float *, float *)
+
+import cython
+from cython.view cimport array as cvarray
+
+import numpy as np
+cimport numpy as np
+
+np.import_array()
+
+ctypedef np.float_t DTYPE_t
+DTYPE = np.float
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+@cython.initializedcheck(False)
+def py_get_Reynolds_SST(np.ndarray[float, ndim=1] lat,
+                        np.ndarray[float, ndim=1] lon,
+                        char *anc_dir, char *sst_file, sst):
+
+    # cdef np.ndarray sst = np.zeros((3232*3200, ), order='C', dtype=np.float32)
+    if not sst.flags['C_CONTIGUOUS']:
+        sst = np.ascontiguousarray(sst)
+
+    cdef float[::1] sst_mv = sst
+
+    get_Reynolds_SST(&lat[0], &lon[0], anc_dir, sst_file, &sst_mv[0])
+
+    return sst
+
+
+def py_get_NDVI_background(np.ndarray[float, ndim=1] lat,
+                           np.ndarray[float, ndim=1] lon,
+                           char *anc_dir, char *ndvi_file, ndvi):
+
+    if not ndvi.flags['C_CONTIGUOUS']:
+        ndvi = np.ascontiguousarray(ndvi)
+
+    cdef float[::1] ndvi_mv = ndvi
+
+    get_NDVI_background(&lat[0], &lon[0], anc_dir, ndvi_file, &ndvi_mv[0])
+
+    return ndvi
+
+
+def py_get_Olson_eco(np.ndarray[float, ndim=1] lat,
+                     np.ndarray[float, ndim=1] lon,
+                     char *anc_dir, eco):
+
+    if not eco.flags['C_CONTIGUOUS']:
+        eco = np.ascontiguousarray(eco)
+
+    cdef unsigned char[::1] eco_mv = eco
+
+    get_Olson_eco(&lat[0], &lon[0], anc_dir, &eco_mv[0])
+
+    return eco
+
+
+def py_get_GEOS(np.ndarray[float, ndim=1] lat, np.ndarray[float, ndim=1] lon, char *startTime,
+                char *anc_dir, char *geos1, char *geos2, char *geos_lnd, char *geos_ocn, char *geos_cnst,
+                tpw, snowfr, icefr, ocnfr, landicefr, sfct):
+
+
+#    for v in geos_data:
+#        if not geos_data[v].flags['C_CONTIGUOUS']:
+#            geos_data[v] = np.ascontiguousarray(geos_data[v])
+    if not tpw.flags['C_CONTIGUOUS']:
+        tpw = np.ascontiguousarray(tpw)
+    if not snowfr.flags['C_CONTIGUOUS']:
+        snowfr = np.ascontiguousarray(snowfr)
+    if not icefr.flags['C_CONTIGUOUS']:
+        icefr = np.ascontiguousarray(icefr)
+    if not ocnfr.flags['C_CONTIGUOUS']:
+        ocnfr = np.ascontiguousarray(ocnfr)
+    if not landicefr.flags['C_CONTIGUOUS']:
+        landicefr = np.ascontiguousarray(landicefr)
+    if not sfct.flags['C_CONTIGUOUS']:
+        sfct = np.ascontiguousarray(sfct)
+
+    cdef float[::1] tpw_mv = tpw  # geos_data['tpw']
+    cdef float[::1] snowfr_mv = snowfr  # geos_data['snowfr']
+    cdef float[::1] icefr_mv = icefr  # geos_data['icefr']
+    cdef float[::1] ocnfr_mv = ocnfr  # geos_data['geos_ocnfr']
+    cdef float[::1] landicefr_mv = landicefr  # geos_data['landicefr']
+    cdef float[::1] sfct_mv = sfct  # geos_data['sfct']
+
+    get_GEOS(&lat[0], &lon[0], startTime, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst,
+             &tpw_mv[0], &snowfr_mv[0], &icefr_mv[0], &ocnfr_mv[0], &landicefr_mv[0], &sfct_mv[0])
+
+#    geos_dict = {'tpw': geos_data['tpw'],
+#                 'snowfr': geos_data['snowfr'],
+#                 'icefr': geos_data['icefr'],
+#                 'geos_ocnfr': geos_data['geos_ocnfr'],
+#                 'landicefr': geos_data['landicefr'],
+#                 'sfct': geos_data['sfct']
+#                 }
+#    geos_dict = {'tpw': tpw,
+#                 'snowfr': snowfr,
+#                 'icefr': icefr,
+#                 'ocnfr': ocnfr,
+#                 'landicefr': landicefr,
+#                 'sfct': sfct
+#                 }
+    print(tpw[:10])
+    return tpw

main.py

 import ruamel_yaml as yml
 import numpy as np
+# import xarray as xr
 
 import read_data as rd
 import tests
 
+import ancillary_data as anc
+
 
 def main():
     datapath = '/ships19/hercules/pveglio/neige_data/snpp_test_input'
@@ -11,8 +14,62 @@ def main():
     fname_geo = 'VNP03MOD.A2014213.1548.001.2017301015705.uwssec.nc'
     thresh_file = '/home/pveglio/mvcm_leo/thresholds/new_thresholds.mvcm.snpp.v1.0.0.yaml'
 
+    anc_dir = '/ships19/hercules/pveglio/neige_data/snpp_test_input/ancillary/2014_08_01_213'
+    sst_file = 'oisst.20140730'
+    ndvi_file = 'NDVI.FM.c004.v2.0.WS.00-04.209.hdf'
+    geos_cnst = 'GEOS.fp.asm.const_2d_asm_Nx.00000000_0000.V01.nc4'
+    geos_lnd = 'GEOS.fpit.asm.tavg1_2d_lnd_Nx.GEOS5124.20140801_1530.V01.nc4'
+    geos_ocn = 'GEOS.fpit.asm.tavg1_2d_ocn_Nx.GEOS5124.20140801_1530.V01.nc4'
+    geos1 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1500.V01.nc4'
+    geos2 = 'GEOS.fpit.asm.inst3_2d_asm_Nx.GEOS5124.20140801_1800.V01.nc4'
+
+    start_time = '2014-08-01 15:48:00.000'
     viirs_data = rd.read_data('viirs', f'{datapath}/{fname_l1b}', f'{datapath}/{fname_geo}')
 
+    sst = np.empty((3232*3200, ), dtype=np.float32)
+    ndvi = np.empty((3232*3200, ), dtype=np.float32)
+    eco = np.empty((3232*3200, ), dtype=np.ubyte)
+
+    geos_data = {'tpw': np.empty((3232*3200, ), dtype=np.float32),
+                 'snowfr': np.empty((3232*3200, ), dtype=np.float32),
+                 'icefr': np.empty((3232*3200, ), dtype=np.float32),
+                 'ocnfr': np.empty((3232*3200, ), dtype=np.float32),
+                 'landicefr': np.empty((3232*3200, ), dtype=np.float32),
+                 'sfct': np.empty((3232*3200, ), dtype=np.float32),
+                 }
+
+    tpw = np.ones((3232*3200, ), dtype=np.float32)
+    snowfr = np.empty((3232*3200, ), dtype=np.float32)
+    icefr = np.empty((3232*3200, ), dtype=np.float32)
+    ocnfr = np.empty((3232*3200, ), dtype=np.float32)
+    landicefr = np.empty((3232*3200, ), dtype=np.float32)
+    sfct = np.empty((3232*3200, ), dtype=np.float32)
+
+    sst = anc.py_get_Reynolds_SST(viirs_data.latitude.values.reshape((3232*3200, )),
+                                  viirs_data.longitude.values.reshape((3232*3200, )),
+                                  anc_dir, sst_file, sst)
+
+    ndvi = anc.py_get_NDVI_background(viirs_data.latitude.values.reshape((3232*3200, )),
+                                      viirs_data.longitude.values.reshape((3232*3200, )),
+                                      anc_dir, ndvi_file, ndvi)
+
+    eco = anc.py_get_Olson_eco(viirs_data.latitude.values.reshape((3232*3200, )),
+                               viirs_data.longitude.values.reshape((3232*3200, )),
+                               anc_dir, eco)
+
+    tpw = anc.py_get_GEOS(viirs_data.latitude.values.reshape((3232*3200, )),
+                          viirs_data.longitude.values.reshape((3232*3200, )),
+                          start_time, anc_dir, geos1, geos2, geos_lnd, geos_ocn, geos_cnst,
+                          tpw, snowfr, icefr, ocnfr, landicefr, sfct)
+
+    print(tpw[:10])
+    #print(geos_data['tpw'][:10])
+    #print(geos_data['snowfr'][:10])
+    #print(geos_data['icefr'][:10])
+    #print(geos_data['ocnfr'][:10])
+    #print(geos_data['landicefr'][:10])
+    #print(geos_data['sfct'][:10])
+
     with open(thresh_file) as f:
         text = f.read()
     thresholds = yml.safe_load(text)