Python SD.end方法代码示例

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def main(cal_file, with_cp):

    from pyhdf.SD import SD

    if with_cp:
        cmd = 'cp %s /home/noel/scratch/' % (cal_file)
        print "running "+cmd
        os.system(cmd)
        filename = os.path.basename(cal_file)
        cal_file = '/home/noel/scratch/' + filename
                        
    print 'Reading ' + cal_file 
    
    vars = ['Latitude', 'Longitude', 
            'Total_Attenuated_Backscatter_532', 'Attenuated_Backscatter_1064', 'Perpendicular_Attenuated_Backscatter_532',
            'Pressure', 'Temperature', 'Molecular_Number_Density', 'Tropopause_Height', 'Surface_Elevation']
    
    hdf = SD(cal_file)
    for var in vars:
        print 'Reading ' + var
        hdf_var = hdf.select(var)
        data = hdf_var.get()
        hdf_var.endaccess()
    hdf.end()
    
    print 'ok.'
    if with_cp:
        print 'Removing '+filename
        cmd = 'rm -f /home/noel/scratch/' + filename
        os.system(cmd)

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def read_rrc(inpath):
    '''Read rrc data m*n from hdf file'''

    '''b1-5;b13-16 for MODIS Rrc
        Rrc_1238 Rrc_443-862 ozone senz solz for VIIRS rrc   
    '''  
    hdf = SD(inpath, SDC.READ)
    #dts = sorted(hdf.datasets().keys())
    modis_key = ['CorrRefl_01','CorrRefl_02','CorrRefl_03','CorrRefl_04','CorrRefl_05',
                 'CorrRefl_13','CorrRefl_14','CorrRefl_15','CorrRefl_16']
    viirs_key = ['Rrc_443','Rrc_486','Rrc_551','Rrc_671','Rrc_745','Rrc_862','Rrc_1238']
    mission = os.path.basename(inpath)[0]
    if mission =='A' or mission =='T':keys = modis_key
    elif mission=='V':keys = viirs_key
    else:keys = hdf.datasets().keys()
    for i,dt in enumerate(keys):
        print(i,dt)
        band = hdf.select(dt)[:,:]        
        if i==0:             
            limit = (band.shape[0],band.shape[1],len(keys))            
            rrc = np.zeros(limit,dtype = np.float)
            rrc[:,:,i] = band
        else:
            rrc[:,:,i] = band
    hdf.end()
    print(rrc.shape)
    return rrc

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def write_interpolated(filename, f0, f1, fact, datasets):
    '''
    interpolate two hdf files f0 and f1 using factor fact, and
    write the result to filename
    '''

    hdf = SD(filename, SDC.WRITE|SDC.CREATE)
    for dataset in datasets:

        try:
            info = SD(f0).select(dataset).info()
        except:
            print >> stderr, 'Error loading %s in %s' % (dataset, f0)
            raise

        typ  = info[3]
        shp  = info[2]
        met0 = SD(f0).select(dataset).get()
        met1 = SD(f1).select(dataset).get()

        interp = (1-fact)*met0 + fact*met1

        interp = interp.astype({
                SDC.INT16: 'int16',
                SDC.FLOAT32: 'float32',
                SDC.FLOAT64: 'float64',
            }[typ])

        # write
        sds = hdf.create(dataset, typ, shp)
        sds[:] = interp[:]
        sds.endaccess()

    hdf.end()

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def read_var_ikslice(filename,var_name,i,k,tc,pc):
    hdffile = SD(filename,SDC.READ)
    ni = hdffile.select('X_grid').ni-1
    nj = hdffile.select('X_grid').nj-1

    if var_name not in ['br','btheta','bphi','vr','vtheta','vphi','bt','bp','vt','vp']:
        var=hdffile.select(var_name+'_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
    else:
        if var_name in ['br','btheta','bphi','bt','bp']:
            bx=hdffile.select('bx_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
            by=hdffile.select('by_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
            bz=hdffile.select('bz_').get(start=(k,0,i),count=(1,nj,1)).squeeze()

            if var_name=='br':
                var     = bx*cos(pc[k])*sin(tc) + by*sin(pc[k])*sin(tc) + bz*cos(tc)
            elif (var_name=='btheta' or var_name=='bt'):
                var = bx*cos(pc[k])*cos(tc) + by*sin(pc[k])*cos(tc) - bz*sin(tc)
            else:
                var   =-bx*sin(pc[k])            + by*cos(pc[k])
        else:
            vx=hdffile.select('vx_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
            vy=hdffile.select('vy_').get(start=(k,0,i),count=(1,nj,1)).squeeze()
            vz=hdffile.select('vz_').get(start=(k,0,i),count=(1,nj,1)).squeeze()

            if var_name=='vr':
                var    = vx*cos(pc[k])*sin(tc) + vy*sin(pc[k])*sin(tc) + vz*cos(tc)
            elif (var_name=='vtheta' or var_name=='vt'):
                var = vx*cos(pc[k])*cos(tc) + vy*sin(pc[k])*cos(tc) - vz*sin(tc)
            else:
                var   =-vx*sin(pc[k])            + vy*cos(pc[k])
    hdffile.end()
    return(var)

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def read_var_islice(filename,var_name,i,thetac,phic):
    nk = phic.size
    nj = thetac.size
    phic = phic[:,None]
    thetac = thetac[None,:]

    hdffile = SD(filename,SDC.READ)
    if var_name not in ['br','btheta','bphi','vr','vtheta','vphi']:
        var=hdffile.select(var_name+'_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
    else:
        if var_name in ['br','btheta','bphi']:
            bx=hdffile.select('bx_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
            by=hdffile.select('by_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
            bz=hdffile.select('bz_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()

            if var_name=='br':
                var     = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
            elif var_name=='btheta':
                var = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
            else:
                var   =-bx*sin(phic)            + by*cos(phic)
        else:
            vx=hdffile.select('vx_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
            vy=hdffile.select('vy_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()
            vz=hdffile.select('vz_').get(start=(0,0,i),count=(nk,nj,1)).squeeze()

            if var_name=='vr':
                var    = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
            elif var_name=='vtheta':
                var = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
            else:
                var   =-vx*sin(phic)            + vy*cos(phic)
    hdffile.end()
    return(var)

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def read_var_point(filename,var_name,i,j,k,thetac,phic):
    thetac = thetac[j]
    phic   = phic[k]


    hdffile = SD(filename,SDC.READ)
    if var_name not in ['br','btheta','bphi','vr','vtheta','vphi']:
        var=hdffile.select(var_name+'_').get(start=(k,j,i),count=(1,1,1)).squeeze()
    else:
#        R,theta,phi=r_theta_phi_uniform(filename)

        if var_name in ['br','btheta','bphi']:
            bx=hdffile.select('bx_').get(start=(k,j,i),count=(1,1,1)).squeeze()
            by=hdffile.select('by_').get(start=(k,j,i),count=(1,1,1)).squeeze()
            bz=hdffile.select('bz_').get(start=(k,j,i),count=(1,1,1)).squeeze()

            if var_name=='br':
                var     = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
            elif var_name=='btheta':
                var = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
            else:
                var   =-bx*sin(phic)            + by*cos(phic)
        else:
            vx=hdffile.select('vx_').get(start=(k,j,i),count=(1,1,1)).squeeze()
            vy=hdffile.select('vy_').get(start=(k,j,i),count=(1,1,1)).squeeze()
            vz=hdffile.select('vz_').get(start=(k,j,i),count=(1,1,1)).squeeze()

            if var_name=='vr':
                var    = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
            elif var_name=='vtheta':
                var = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
            else:
                var   =-vx*sin(phic)            + vy*cos(phic)
    hdffile.end()
    return(var)

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def readMOD35L2(fname, geoloc_only=False):
    hdf_file = SD(HDFDIR + fname)
    if not geoloc_only:
        cloud_mask = hdf_file.select('Cloud_Mask').get()
    lon = hdf_file.select('Longitude').get()
    lat = hdf_file.select('Latitude').get()
    hdf_file.end()
    
    if not geoloc_only:
        cld_msk = uint8(cloud_mask[0])
        cloud = cld_msk & 6 # 0, 2, 4, 6
        land = cld_msk & 192 # 0, 64, 128, 192
    
        cloud[cloud==0] = 1 # 0 -> Cloud
        cloud[cloud!=1] = 0 # 2, 4, 6 -> No cloud

        coast = land
        coast[coast==64] = 1 # 64 -> Coast
        coast[coast!=1] = 0 # 0, 128, 192 -> Not coast

        land[land!=0] = 1 # 64, 128, 192 -> Land, 0 -> Water
        
        return lon, lat, cloud, land, coast

    return lon, lat

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
    def read(self, filename, **kwargs):
        """Read the data"""
        from pyhdf.SD import SD
        import datetime

        #print "*** >>> Read the hdf-eos file!"
        root = SD(filename)
    
        # Get all the Attributes:
        # Common Attributes, Data Time,
        # Data Structure and Scene Coordinates
        for key in root.attributes().keys():
            self._eoshdf_info[key] = root.attributes()[key]

        # Start Time - datetime object
        starttime = datetime.datetime.strptime(self._eoshdf_info['Start Time'][0:13], 
                                               "%Y%j%H%M%S")
        msec = float(self._eoshdf_info['Start Time'][13:16])/1000.
        self.starttime = starttime + datetime.timedelta(seconds=msec)
    
        # End Time - datetime object
        endtime = datetime.datetime.strptime(self._eoshdf_info['End Time'][0:13], 
                                             "%Y%j%H%M%S")
        msec = float(self._eoshdf_info['End Time'][13:16])/1000.
        self.endtime = endtime + datetime.timedelta(seconds=msec)

        # What is the leading 'H' doing here?
        sensor_name = self._eoshdf_info['Sensor Name'][1:-1].lower()
        try:
            self.satid = EOS_SATELLITE[sensor_name]
        except KeyError:
            LOG.error("Failed setting the satellite id - sat-name = ", 
                      sensor_name)
            
        self.orbit = self._eoshdf_info['Orbit Number']
        self.shape = (self._eoshdf_info['Number of Scan Control Points'],
                      self._eoshdf_info['Number of Pixel Control Points'])

        #try:
        if 1:
            value = root.select(self.name)
            attr = value.attributes()
            data = value.get()

            self.attr = attr
            band = data
            if self.name in FLAGS_QUALITY:
                self.data = band
            else:
                nodata = attr['bad_value_scaled']
                self.data = (np.ma.masked_equal(band, nodata) * 
                             attr['slope'] + attr['intercept'])
            
            value.endaccess()
        #except:
        #    pass

        root.end()
        self.filled= True

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def readhdf(filename, fieldname, ignoreE = True):
    try:
        hdf = SD(filename, SDC.READ)
        data = hdf.select(fieldname)[:].copy()
        hdf.end()
    except Exception, e:
        if not ignoreE:print e
        data = Dataset(filename)[fieldname][:].copy()

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def read_var(fname,varname,normalized=False):
    f     = SD(fname,SDC.READ)
    phi   = f.select('fakeDim0')[:]
    theta = f.select('fakeDim1')[:]
    r     = f.select('fakeDim2')[:]
    var   = f.select('Data-Set-2')[:]
    f.end()

    if normalized:
        return(phi,theta,r,var)
    else:
        return(phi,theta,r*mas_units['length'],var*mas_units[varname])

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def t12(mfile):
    hdffile = SD(mfile)
    x = hdffile.select('MYD021KM_EV_1KM_Emissive_Band32')
    scale_factor = x.attributes()['scale_factor']
    add_offset = x.attributes()['add_offset']
    x = x[:]
    hdffile.end()    
    
    idx = x > 65534
    x = (x - add_offset) * scale_factor
    t = planck (12.02, x)
    t[idx] = np.nan
    return t

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
class GeoProf(object):
    
    def __init__(self, filename):
        
        self.hdf = SD(filename)
        self.filename = filename
        self.z = filename[-2:]
        self.orbit = filename[-15:-4]
        self.id = filename[-25:-4]
        
    def close(self):
        self.hdf.end()
        self.hdf = None
        
    def _read_var(self, varname):
        hdfvar = self.hdf.select(varname)
        data = hdfvar[:]
        hdfvar.endaccess()
        return data
    
    def coords(self):
        lat = self._read_var('Latitude')
        lon = self._read_var('Longitude')
        return lon, lat
    
    def time(self):
        time = self._read_var('Time')
        return time
    
    
    def altitude(self):
        ''' altitude in kilometers '''
        
        alt = self._read_var('Height') / 1e3
        return alt
    
    def cloudmask(self):
        '''
        "0 = No cloud detected\n",
        "1 = likely bad data\n",
        "5 = likely ground clutter\n",
        "5-10 = week detection found using along track integration\n",
        "20 to 40 = Cloud detected .. increasing values represents clouds with lower chance of a being a false detection" ;
        _FillValue = '\200' ;
        
        shape [nprof, nalt]
        '''
        
        cm = self._read_var('CPR_Cloud_mask')
        
        return cm

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def read(filename):
    """
    Returns R,theta,phi,Rc,thetac,phic,br,btheta,bphi,vr,vtheta,vphi,rho,cs
    """

    hdffile = SD(filename,SDC.READ)

    x=hdffile.select('X_grid').get()
    y=hdffile.select('Y_grid').get()
    z=hdffile.select('Z_grid').get()
    bx=hdffile.select('bx_').get()[:-1,:-1,:-1]
    by=hdffile.select('by_').get()[:-1,:-1,:-1]
    bz=hdffile.select('bz_').get()[:-1,:-1,:-1]
    vx=hdffile.select('vx_').get()[:-1,:-1,:-1]
    vy=hdffile.select('vy_').get()[:-1,:-1,:-1]
    vz=hdffile.select('vz_').get()[:-1,:-1,:-1]
    rho=hdffile.select('rho_').get()[:-1,:-1,:-1]
    cs=hdffile.select('c_').get()[:-1,:-1,:-1]

    t=hdffile.time
    hdffile.end()

    # =========== Cell centers ==============
    xc=0.125*( x[:-1,:-1,:-1]+x[1:,:-1,:-1]+x[:-1,1:,:-1]+x[:-1,:-1,1:]+
              x[1:,1:,:-1]+x[1:,:-1,1:]+x[:-1,1:,1:]+x[1:,1:,1:] )
    yc=0.125*( y[:-1,:-1,:-1]+y[1:,:-1,:-1]+y[:-1,1:,:-1]+y[:-1,:-1,1:]+
              y[1:,1:,:-1]+y[1:,:-1,1:]+y[:-1,1:,1:]+y[1:,1:,1:] )
    zc=0.125*( z[:-1,:-1,:-1]+z[1:,:-1,:-1]+z[:-1,1:,:-1]+z[:-1,:-1,1:]+
              z[1:,1:,:-1]+z[1:,:-1,1:]+z[:-1,1:,1:]+z[1:,1:,1:] )
    # =======================================
    
    R=sqrt(x**2+y**2+z**2)
    theta=arccos(z/R)
    phi=arctan2(y,x)
    phi[phi<0]+=2*pi

    Rc=sqrt(xc**2+yc**2+zc**2)
    thetac=arccos(zc/Rc)
    phic=arctan2(yc,xc)
    phic[phic<0]+=2*pi


    br     = bx*cos(phic)*sin(thetac) + by*sin(phic)*sin(thetac) + bz*cos(thetac)
    btheta = bx*cos(phic)*cos(thetac) + by*sin(phic)*cos(thetac) - bz*sin(thetac)
    bphi   =-bx*sin(phic)            + by*cos(phic)

    vr     = vx*cos(phic)*sin(thetac) + vy*sin(phic)*sin(thetac) + vz*cos(thetac)
    vtheta = vx*cos(phic)*cos(thetac) + vy*sin(phic)*cos(thetac) - vz*sin(thetac)
    vphi   =-vx*sin(phic)            + vy*cos(phic)

    return(t,R,theta,phi,Rc,thetac,phic,br,btheta,bphi,vr,vtheta,vphi,rho,cs)

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
def callback(body, message):
    """Do actual work."""

    logger.info("body in callback() is %s" % body)

    # pull lat/lon, time
    path = body
    sd = SD(path)
    lat = N.array(sd.select('Latitude').get())
    lon = N.array(sd.select('Longitude').get())
    t = N.array(sd.select('Time').get())
    sd.end()
    #logger.info("lat: %s" % str(lat.shape))
    #logger.info("lon: %s" % str(lon.shape))
    #logger.info("time: %s" % str(t.shape))

    # build metadata json
    id = os.path.basename(path)
    md = {
        "id": id,
        "dataset": "AIRX2RET",
        "starttime": t[0,0],
        "endtime": t[44,29],
        "location": {
            "coordinates": [[
                [ lon[0,0], lat[0,0] ],
                [ lon[0,29], lat[0,29] ],
                [ lon[44,29], lat[44,29] ],
                [ lon[44,0], lat[44,0] ],
                [ lon[0,0], lat[0,0] ],
            ]], 
            "type": "polygon"
        }, 
        "urls": "http://mozart/data/public/products/%s" % id
    }

    # publish
    pub_dir = '/data/public/products'
    ensure_dir(pub_dir)
    shutil.move(path, os.path.join(pub_dir, id))

    # insert into ElasticSearch
    index = doctype = 'airs'
    conn = ES('http://localhost:9200')
    mapping = json.load(open('grq_mapping.json'))
    if not conn.indices.exists_index(index):
        conn.indices.create_index(index, mapping)
    conn.indices.put_mapping(doctype, mapping, index)
    ret = conn.index(md, index, doctype, md['id'])

    message.ack()

# 需要导入模块: from pyhdf.SD import SD [as 别名]
# 或者: from pyhdf.SD.SD import end [as 别名]
    def parseMetadata(self, filepath):

        metadata = {}
    
        dir, filename = os.path.split(filepath)
        if re.match(FILENAME_PATTERN, filename):
            logging.info("Parsing HDF file=%s" % filepath)

            # open HDF file
            try:
                hdfFile = SD(filepath, SDC.READ)
            except HDF4Error as e:
                logging.info(e)
                raise e

            # variables
            variables = hdfFile.datasets().keys()

            # time fields
            year = hdfFile.select('Year')[:]
            month = hdfFile.select('Month')[:]
            day = hdfFile.select('Day')[:]
            hour = hdfFile.select('Hour')[:]
            minute = hdfFile.select('Minute')[:]
            second = hdfFile.select('Seconds')[:]

            # space fields
            lon = hdfFile.select('Longitude')[:]
            lat = hdfFile.select('Latitude')[:]

            datetimes = []
            lats = []
            lons = []
            for t in range(22):
                for x in range(15):
                    if year[t,x] != -9999:

                        datetimes.append( dt.datetime(year[t,x],month[t,x],day[t,x],hour[t,x],minute[t,x],second[t,x], tzinfo=tzutc()) )
                        lons.append( lon[t,x] )
                        lats.append( lat[t,x] )
                        
            # store metadata values
            storeMetadata(metadata, np.asarray(lons), np.asarray(lats), np.asarray(datetimes), variables)

            # close HDF file
            hdfFile.end()

        return metadata