3、利用python读取Landsat系列数据的元数据(*MTL.txt文件) 在前几年NASA里的某团对发布了处理Landsat的python(2.7)工具包（ dnppy ），但我在电脑上安装了多次都没有成功，于是放弃安装了。根据其官方介绍，该工具包能读元数据(*MTL.txt文件)、能实现反射率定标(包括星上(TOA)反射率和地表(surface)反射率)和辐射定标，感觉很实用。我又不愿意自己去写，于是就在他们给的代码了修改了一下，以为己用。这里直接贴上了。将该代码复制保存为.py文件(如lst_meta.py)，然后再import就可以直接调用里面的函数了(import lst_meta) __author__ = 'jwely(modified by wenzhaofei)' # standard imports from datetime import datetime, timedelta from numpy import radians, ndarray, sin, cos, degrees, arctan2, arcsin, tan, arccos import inspect class solar: Object class for handling solar calculations. Many equations are taken from the excel sheet at this url : It requires a physical location on the earth and a datetime object :param lat: decimal degrees latitude (float OR numpy array) :param lon: decimal degrees longitude (float OR numpy array) :param time_zone: float of time shift from GMT (such as -5 for EST) :param date_time_obj: either a timestamp string following fmt or a datetime obj :param fmt: if date_time_obj is a string, fmt is required to interpret it :param slope: slope of land at lat,lon for solar energy calculations :param aspect: aspect of land at lat,lon for solar energy calculations An instance of this class may have the following attributes: =================== =========================================== ======== attribute description type =================== =========================================== ======== lat latitude (array) lon longitude (array) tz time zone (scalar) rdt reference datetime object (date_time_obj) (scalar) slope slope, derivative of DEM (array) aspect aspect (north is 0, south is 180) (array) ajd absolute julian day (scalar) ajc absolute julian century (scalar) geomean_long geometric mean longitude of the sun (scalar) geomean_anom geometric mean longitude anomaly of the sun (scalar) earth_eccent eccentricity of earths orbit (scalar) sun_eq_of_center the suns equation of center (scalar) true_long true longitude of the sun (scalar) true_anom true longitude anomaly of the sun (scalar) app_long the suns apparent longitude (scalar) oblique_mean_elip earth oblique mean ellipse (scalar) oblique_corr correction to earths oblique elipse (scalar) right_ascension suns right ascension angle (scalar) declination solar declination angle (scalar) equation_of_time equation of time (minutes) (scalar) hour_angle_sunrise the hour angle at sunrise (array) solar_noon LST of solar noon (array) sunrise LST of sunrise time (array) sunset LST of sunset time (array) sunlight LST fractional days of sunlight (array) true_solar LST for true solar time (array) hour_angle total hour angle (array) zenith zenith angle (array) elevation elevation angle (array) azimuth azimuthal angle (array) rad_vector radiation vector (distance in AU) (scalar) earth_distance earths distance to sun in meters (scalar) norm_irradiance incident solar energy at earth distance (scalar) =================== =========================================== ======== Units used by this class unless otherwise labeled - angle = degrees - distance = meters - energy = watts or joules - time = mostly in datetime objects. labeled in most cases. Planned improvements 1. DONE. Inputs of numpy arrays for lat and lon needs to be allowed. 2. inputs of a numpy array DEM for slope/aspect effects on incident solar energy Present performance To process about one landsat tile (7300^2 matrix) requires 9GB of memory and takes 45 seconds to process on a single 3.3GHz thread. It would be nice to get the same output to run on ~5GB of memory so a 8GB system could handle it. def __init__(self, lat, lon, date_time_obj, time_zone = 0, fmt = False, slope = None, aspect = None): Initializes critical spatial and temporal information for solar object. # empty list of class attributes self.ajc = None # abs julian century (defined on __init__) self.ajd = None # abs julian day (defined on __init__) self.app_long = None self.atmo_refraction = None self.azimuth = None self.declination = None self.earth_distance = None self.earth_eccent = None self.elevation = None self.elevation_noatmo = None self.equation_of_time = None self.frac_day = None self.geomean_anom = None self.geomean_long = None self.hour_angle = None self.hour_angle_sunrise = None self.lat = lat # lattitude (E positive)- float self.lat_r = radians(lat) # lattitude in radians self.lon = lon # longitude (N positive)- float self.lon_r = radians(lon) # longitude in radians self.norm_irradiance = None self.oblique_corr = None self.oblique_mean_elip = None self.rad_vector = None self.rdt = None # reference datetime (defined on __init__) self.right_ascension = None self.solar_noon = None self.solar_noon_time = None self.sun_eq_of_center = None self.sunlight = None self.sunlight_time = None self.sunrise = None self.sunrise_time = None self.sunset = None self.sunset_time = None self.true_anom = None self.true_long = None self.true_solar = None self.true_solar_time = None self.tz = None # time zone (defined on __init__) self.zenith = None # slope and aspect self.slope = slope self.aspect = aspect # Constants as attributes self.sun_surf_rad = 63156942.6 # radiation at suns surface (W/m^2) self.sun_radius = 695800000. # radius of the sun in meters self.orbital_period = 365.2563630 # num of days it takes earth to revolve self.altitude = -0.01448623 # altitude of center of solar disk # sets up the object with some subfunctions self._set_datetime(date_time_obj, fmt, GMT_hour_offset = time_zone) # specify if attributes are scalar floats or numpy array floats if isinstance(lat, ndarray) and isinstance(lon, ndarray): self.is_numpy = True else: self.is_numpy = False return def _set_datetime(self, date_time_obj, fmt = False, GMT_hour_offset = 0): sets the critical time information including absolute julian day/century. Accepts datetime objects or a datetime string with format :param date_time_obj: datetime object for time of solar calculations. Will also accept string input with matching value for fmt param :param fmt: if date_time_obj is input as a string, fmt allows it to be interpreted :param GMT_hour_offset: Number of hours from GMT for timezone of calculation area. # if input is datetime_obj set it if isinstance(date_time_obj, datetime): self.rdt = date_time_obj self.rdt += timedelta(hours = -GMT_hour_offset) elif isinstance(date_time_obj, str) and isinstance(fmt, str): self.rdt = datetime.strptime(date_time_obj,fmt) self.rdt += timedelta(hours = -GMT_hour_offset) else: raise Exception(bad datetime!) self.tz = GMT_hour_offset # uses the reference day of january 1st 2000 jan_1st_2000_jd = 2451545 jan_1st_2000 = datetime(2000,1,1,12,0,0) time_del = self.rdt - jan_1st_2000 self.ajd = float(jan_1st_2000_jd) + float(time_del.total_seconds())/86400 self.ajc = (self.ajd - 2451545)/36525.0 return def get_geomean_long(self): :return geomean_long: geometric mean longitude of the sun if not self.geomean_long is None: return self.geomean_long self.geomean_long = (280.46646 + self.ajc * (36000.76983 + self.ajc*0.0003032)) % 360 return self.geomean_long def get_geomean_anom(self): calculates the geometric mean anomoly of the sun if not self.geomean_anom is None: return self.geomean_anom self.geomean_anom = (357.52911 + self.ajc * (35999.05029 - 0.0001537 * self.ajc)) return self.geomean_anom def get_earth_eccent(self): :return earth_eccent: precise eccentricity of earths orbit at referece datetime if not self.earth_eccent is None: return self.earth_eccent self.earth_eccent = 0.016708634 - self.ajc * (4.2037e-5 + 1.267e-7 * self.ajc) return self.earth_eccent def get_sun_eq_of_center(self): :return sun_eq_of_center: the suns equation of center if not self.sun_eq_of_center is None: return self.sun_eq_of_center if self.geomean_anom is None: self.get_geomean_anom() ajc = self.ajc gma = radians(self.geomean_anom) self.sun_eq_of_center = sin(gma) * (1.914602 - ajc*(0.004817 + 0.000014 * ajc)) + \\ sin(2*gma) * (0.019993 - 0.000101 * ajc) + \\ sin(3*gma) * 0.000289 return self.sun_eq_of_center def get_true_long(self): :return true_long: the tru longitude of the sun if not self.true_long is None: return self.true_long if self.geomean_long is None: self.get_geomean_long() if self.sun_eq_of_center is None: self.get_sun_eq_of_center() self.true_long = self.geomean_long + self.sun_eq_of_center return self.true_long def get_true_anom(self): :return true_anom: calculates the true anomaly of the sun if not self.true_anom is None: return self.true_anom if self.geomean_long is None: self.get_geomean_long() if self.sun_eq_of_center is None: self.get_sun_eq_of_center() self.true_anom = self.geomean_anom + self.sun_eq_of_center return self.true_anom def get_rad_vector(self): :return rad_vector: incident radiation vector to surface at ref_datetime (AUs) if not self.rad_vector is None: return self.rad_vector if self.earth_eccent is None: self.get_earth_eccent() if self.true_anom is None: self.get_true_anom() ec = self.earth_eccent ta = radians(self.true_anom) self.rad_vector = (1.000001018*(1 - ec**2)) / (1 + ec *cos(ta)) return self.rad_vector def get_app_long(self): :return app_long: calculates apparent longitude of the sun if not self.app_long is None: return self.app_long if self.true_long is None: self.get_true_long() stl = self.true_long ajc = self.ajc self.app_long = stl - 0.00569 - 0.00478 * sin(radians(125.04 - 1934.136 * ajc)) return self.app_long def get_oblique_mean_elip(self): :return oblique_mean_elip: oblique mean elliptic of earth orbit if not self.oblique_mean_elip is None: return self.oblique_mean_elip ajc = self.ajc self.oblique_mean_elip = 23 + (26 + (21.448 - ajc * (46.815 + ajc * (0.00059 - ajc * 0.001813)))/60)/60 return self.oblique_mean_elip def get_oblique_corr(self): :return oblique_corr: the oblique correction if not self.oblique_corr is None: return self.oblique_corr if self.oblique_mean_elip is None: self.get_oblique_mean_elip() ome = self.oblique_mean_elip ajc = self.ajc self.oblique_corr = ome + 0.00256 * cos(radians(125.04 - 1934.136 * ajc)) return self.oblique_corr def get_right_ascension(self): :return right_ascension: the suns right ascension angle if not self.right_ascension is None: return self.right_ascension if self.app_long is None: self.get_app_long() if self.oblique_corr is None: self.get_oblique_corr() sal = radians(self.app_long) oc = radians(self.oblique_corr) self.right_ascension = degrees(arctan2(cos(oc) * sin(sal), cos(sal))) return self.right_ascension def get_declination(self): :return declination: solar declination angle at ref_datetime if not self.declination is None: return self.declination if self.app_long is None: self.get_app_long() if self.oblique_corr is None: self.get_oblique_corr() sal = radians(self.app_long) oc = radians(self.oblique_corr) self.declination = degrees(arcsin((sin(oc) * sin(sal)))) return self.declination def get_equation_of_time(self): :return equation_of_time: the equation of time in minutes if not self.equation_of_time is None: return self.equation_of_time if self.oblique_corr is None: self.get_oblique_corr() if self.geomean_long is None: self.get_geomean_long() if self.geomean_anom is None: self.get_geomean_anom() if self.earth_eccent is None: self.get_earth_eccent() oc = radians(self.oblique_corr) gml = radians(self.geomean_long) gma = radians(self.geomean_anom) ec = self.earth_eccent vary = tan(oc/2)**2 self.equation_of_time = 4 * degrees(vary * sin(2*gml) - 2 * ec * sin(gma) + 4 * ec * vary * sin(gma) * cos(2 * gml) - 0.5 * vary * vary * sin(4 * gml) - 1.25 * ec * ec * sin(2 * gma)) return self.equation_of_time def get_hour_angle_sunrise(self): :return hour_angle_sunrise: the hour angle of sunrise if not self.hour_angle_sunrise is None: return self.hour_angle_sunrise if self.declination is None: self.get_declination() d = radians(self.declination) lat = self.lat_r self.hour_angle_sunrise = degrees(arccos((cos(radians(90.833)) / (cos(lat) * cos(d)) - tan(lat) * tan(d)))) return self.hour_angle_sunrise def get_solar_noon(self): :return solar_noon: solar noon in (local sidereal time LST) if not self.solar_noon is None: return self.solar_noon if self.equation_of_time is None: self.get_equation_of_time() lon = self.lon eot = self.equation_of_time tz = self.tz self.solar_noon = (720 - 4 * lon - eot + tz * 60)/1440 # format this as a time for display purposes (Hours:Minutes:Seconds) if self.is_numpy: self.solar_noon_time = timedelta(days = self.solar_noon.mean()) else: self.solar_noon_time = timedelta(days = self.solar_noon) return self.solar_noon def get_sunrise(self): :return sunrise: returns the time of sunrise if not self.sunrise is None: return self.sunrise if self.solar_noon is None: self.get_solar_noon() if self.hour_angle_sunrise is None: self.get_hour_angle_sunrise() sn = self.solar_noon ha = self.hour_angle_sunrise self.sunrise = (sn * 1440 - ha * 4)/1440 # format this as a time for display purposes (Hours:Minutes:Seconds) if self.is_numpy: self.sunrise_time = timedelta(days = self.sunrise.mean()) else: self.sunrise_time = timedelta(days = self.sunrise) return self.sunrise def get_sunset(self): :return sunset: returns the time of sunset if not self.sunset is None: return self.sunset if self.solar_noon is None: self.get_solar_noon() if self.hour_angle_sunrise is None: self.get_hour_angle_sunrise() sn = self.solar_noon ha = self.hour_angle_sunrise self.sunset = (sn * 1440 + ha * 4)/1440 # format this as a time for display purposes (Hours:Minutes:Seconds) if self.is_numpy: self.sunset_time = timedelta(days = self.sunset.mean()) else: self.sunset_time = timedelta(days = self.sunset) return self.sunset def get_sunlight(self): :return sunlight: amount of daily sunlight in fractional days if not self.sunlight is None: return self.sunlight if self.hour_angle_sunrise is None: self.get_hour_angle_sunrise() self.sunlight = 8 * self.hour_angle_sunrise / (60 * 24) # format this as a time for display purposes (Hours:Minutes:Seconds) if self.is_numpy: self.sunlight_time = timedelta(days = self.sunlight.mean()) else: self.sunlight_time = timedelta(days = self.sunlight) return self.sunlight def get_true_solar(self): :return true_solar: true solar time at ref_datetime if not self.true_solar is None: return self.true_solar if self.equation_of_time is None: self.get_equation_of_time lon = self.lon eot = self.equation_of_time # turn reference datetime into fractional days frac_sec = (self.rdt - datetime(self.rdt.year, self.rdt.month, self.rdt.day)).total_seconds() frac_hr = frac_sec / (60 * 60) + self.tz frac_day = frac_hr / 24 self.frac_day = frac_day # now get true solar time self.true_solar = (frac_day * 1440 + eot + 4 * lon - 60 * self.tz) % 1440 # format this as a time for display purposes (Hours:Minutes:Seconds) if self.is_numpy: self.true_solar_time = timedelta(days = self.true_solar.mean() / (60*24)) else: self.true_solar_time = timedelta(days = self.true_solar / (60*24)) return self.true_solar def get_hour_angle(self): :return hour_angle: returns hour angle at ref_datetime if not self.hour_angle is None: return self.hour_angle if self.true_solar is None: self.get_true_solar() ts = self.true_solar # matrix hour_angle calculations if self.is_numpy: ha = ts ha = ha /4 + 180 ha = ha /4 - 180 self.hour_angle = ha # scalar hour_angle calculations else: if ts = 0: self.hour_angle = ts/4 + 180 else: self.hour_angle = ts/4 - 180 return self.hour_angle def get_zenith(self): :return zenith: returns solar zenith angle at ref_datetime if not self.zenith is None: return self.zenith if self.declination is None: self.get_declination() if self.hour_angle is None: self.get_hour_angle() d = radians(self.declination) ha = radians(self.hour_angle) lat = self.lat_r self.zenith = degrees(arccos(sin(lat) * sin(d) + cos(lat) * cos(d) * cos(ha))) return self.zenith def get_elevation(self): :return elevation: returns solar elevation angle at ref_datetime if not self.elevation is None: return self.elevation if self.zenith is None: self.get_zenith() # perform an approximate atmospheric refraction correction # matrix hour_angle calculations # these equations are hideous, but im not sure how to improve them without # adding computational complexity if self.is_numpy: e = 90 - self.zenith ar = e * 0 ar = 0 ar = 58.1 / tan(radians(e )) - \\ 0.07 / tan(radians(e ))**3 + \\ 0.000086 / tan(radians(e ))**5 ar = 1735 + e * \\ (103.4 + e * (-12.79 + e * 0.711)) ar = -20.772 / tan(radians(e )) # scalar hour_angle calculations else: e = 90 - self.zenith er = radians(e) if e 85: ar = 0 elif e 5: ar = 58.1 / tan(er) - 0.07 / tan(er)**3 + 0.000086 / tan(er)**5 elif e -0.575: ar = 1735 + e * (103.4 + e * ( -12.79 + e * 0.711)) else: ar = -20.772 / tan(er) self.elevation_noatmo = e self.atmo_refraction = ar / 3600 self.elevation = self.elevation_noatmo + self.atmo_refraction return self.elevation def get_azimuth(self): :return azimuth: returns solar azimuth angle at ref_datetime if not self.azimuth is None: return self.azimuth if self.declination is None: self.get_declination() if self.hour_angle is None: self.get_hour_angle() if self.zenith is None: self.get_zenith() lat = self.lat_r d = radians(self.declination) ha = radians(self.hour_angle) z = radians(self.zenith) # matrix azimuth calculations # these equations are hideous monsters, but im not sure how to improve them without # adding computational complexity. if self.is_numpy: az = ha * 0 az = (degrees(arccos(((sin(lat ) * cos(z )) - sin(d)) / (cos(lat ) * sin(z )))) + 180) % 360 az = (540 - degrees(arccos(((sin(lat ) * cos(z )) -sin(d))/ (cos(lat ) * sin(z ))))) % 360 self.azimuth = az else: if ha 0: self.azimuth = (degrees(arccos(((sin(lat) * cos(z)) - sin(d)) / (cos(lat) * sin(z)))) + 180) % 360 else: self.azimuth = (540 - degrees(arccos(((sin(lat) * cos(z)) -sin(d))/ (cos(lat) * sin(z))))) % 360 return self.azimuth def get_earth_distance(self): :return earth_distance: distance between the earth and the sun at ref_datetime if self.rad_vector is None: self.get_rad_vector() # convert rad_vector length from AU to meters self.earth_distance = self.rad_vector * 149597870700 return self.earth_distance def get_norm_irradiance(self): Calculates incoming solar energy in W/m^2 to a surface normal to the sun. inst_irradiance is calculated as = sun_surf_radiance\\*(sun_radius / earth_distance)^2 and is then corrected as a function of solar incidence angle :return norm_irradiance: the normal irradiance in W/m^2 if not self.norm_irradiance is None: return self.norm_irradiance if self.earth_distance is None: self.get_earth_distance() ed = self.earth_distance # calculate irradiance to normal surface at earth distance self.norm_irradiance = self.sun_surf_rad * (self.sun_radius / ed)**2 return self.norm_irradiance def get_inc_irradiance(self): calculates the actual incident solar irradiance at a given lat,lon coordinate with adjustments for slope and aspect if they have been given. Not finished. print(this function is unfinished!) return def summarize(self): prints attribute list and corresponding values for key in sorted(self.__dict__.iterkeys()): print({0} {1}.format(key.ljust(20),sc.__dict__ )) return def compute_all(self): Computes and prints all the attributes of this solar object. Spatial averages are printed for numpy array type attributes. print(=*50) print(Interogation of entire matrix of points.) print(Some values displayed below are spatial averages) print(=*50) if self.is_numpy: # print means of lat/lon arrays print(latitude, longitude \\t{0}, {1}.format(self.lat.mean(), self.lon.mean())) else: print(latitude, longitude \\t{0}, {1}.format(self.lat, self.lon)) print(datetime \\t\\t{0} (GMT).format(self.rdt)) print(time zone \\t\\t{0} (GMT offset).format(self.tz)) print() print(abs julian day \\t\\t{0}\\t (day).format(self.ajd)) print(abs julian century \\t{0}\\t (cen).format(self.ajc)) print(suns geomean long \\t{0}\\t (deg).format(self.get_geomean_long())) print(suns geomean anom \\t{0}\\t (deg).format(self.get_geomean_anom())) print(earth eccentricity \\t{0}.format(self.get_earth_eccent())) print(suns eq of center \\t{0}.format(self.get_sun_eq_of_center())) print(suns true long \\t\\t{0}\\t (deg).format(self.get_true_long())) print(suns true anom \\t\\t{0}\\t (deg).format(self.get_true_anom())) print(suns apparent long \\t{0}\\t (deg).format(self.get_app_long())) print(earth obliq mean elip \\t{0}\\t (deg).format(self.get_oblique_mean_elip())) print(earth obliq correction\\t{0}\\t (deg).format(self.get_oblique_corr())) print(sun right ascension \\t{0}\\t (deg).format(self.get_right_ascension())) print(solar declination angle {0}\\t (deg).format(self.get_declination())) print(equation of time \\t{0}\\t (min).format(self.get_equation_of_time)) if self.is_numpy: # print means of hour angle array print(hour angle sunrise\\t{0}\\t (deg).format(self.get_hour_angle_sunrise().mean())) else: print(hour angle sunrise\\t{0}\\t (deg).format(self.get_hour_angle_sunrise())) print() self.get_solar_noon() print(solar noon \\t\\t{0}\\t (HMS - LST).format(self.solar_noon_time)) self.get_sunrise() print(sunrise \\t\\t{0}\\t (HMS - LST).format(self.sunrise_time)) self.get_sunset() print(sunset \\t\\t{0}\\t (HMS - LST).format(self.sunset_time)) self.get_sunlight() print(sunlight durration \\t{0}\\t (HMS).format(self.sunlight_time)) self.get_true_solar() print(true solar time \\t{0}\\t (HMS - LST).format(self.true_solar_time)) print() if self.is_numpy: # print means of these array objects print(hour angle \\t\\t{0}\\t (deg).format(self.get_hour_angle().mean())) print(solar zenith angle \\t{0}\\t (deg).format(self.get_zenith().mean())) print(solar elevation angle \\t{0}\\t (deg).format(self.get_elevation().mean())) print(solar azimuth angle \\t{0}\\t (deg).format(self.get_azimuth().mean())) else: print(hour angle \\t\\t{0}\\t (deg).format(self.get_hour_angle())) print(solar zenith angle \\t{0}\\t (deg).format(self.get_zenith())) print(solar elevation angle \\t{0}\\t (deg).format(self.get_elevation())) print(solar azimuth angle \\t{0}\\t (deg).format(self.get_azimuth())) print() print(radiation vector \\t{0}\\t (AU).format(self.get_rad_vector())) print(earth sun distance \\t{0}(m).format(self.get_earth_distance())) print(norm irradiance \\t{0}\\t (W/m*m).format(self.get_norm_irradiance())) print(=*50) # testing if __name__ == __main__: # use the current time and my time zone my_datestamp = 20150515-120000 # date stamp my_fmt = %Y%m%d-%H%M%S # datestamp format my_tz = -4 # timezone (GMT/UTC) offset my_lat = 37 # lat (N positive) my_lon = -76.4 # lon (E positive) sc = solar(my_lat, my_lon, my_datestamp, my_tz, my_fmt) sc.compute_all() sc.summarize() class landsat_metadata: A landsat metadata object. This class builds is attributes from the names of each tag in the xml formatted .MTL files that come with landsat data. So, any tag that appears in the MTL file will populate as an attribute of landsat_metadata. You can access explore these attributes by using, for example .. code-block:: python from dnppy import landsat meta = landsat.landsat_metadata(my_filepath) # create object from pprint import pprint # import pprint pprint(vars(m)) # pretty print output scene_id = meta.LANDSAT_SCENE_ID # access specific attribute :param filename: the filepath to an MTL file. def __init__(self, filename): There are several critical attributes that keep a common naming convention between all landsat versions, so they are initialized in this class for good record keeping and reference # custom attribute additions self.FILEPATH = filename self.DATETIME_OBJ = None # product metadata attributes self.LANDSAT_SCENE_ID = None self.DATA_TYPE = None self.ELEVATION_SOURCE = None self.OUTPUT_FORMAT = None self.SPACECRAFT_ID = None self.SENSOR_ID = None self.WRS_PATH = None self.WRS_ROW = None self.NADIR_OFFNADIR = None self.TARGET_WRS_PATH = None self.TARGET_WRS_ROW = None self.DATE_ACQUIRED = None self.SCENE_CENTER_TIME = None # image attributes self.CLOUD_COVER = None self.IMAGE_QUALITY_OLI = None self.IMAGE_QUALITY_TIRS = None self.ROLL_ANGLE = None self.SUN_AZIMUTH = None self.SUN_ELEVATION = None self.EARTH_SUN_DISTANCE = None # calculated for Landsats before 8. # read the file and populate the MTL attributes self._read(filename) def _read(self, filename): reads the contents of an MTL file # if the filename input is actually already a metadata class object, return it back. if inspect.isclass(filename): return filename fields = metafile = open(filename, 'r') metadata = metafile.readlines() for line in metadata: # skips lines that contain bad flags denoting useless data AND lines # greater than 1000 characters. 1000 character limit works around an odd LC5 # issue where the metadata has 40,000+ characters of whitespace bad_flags = if not any(x in line for x in bad_flags) and len(line) = 1000: try: line = line.replace( , ) line = line.replace(\\n, ) field_name, field_value = line.split(' = ') fields.append(field_name) values.append(field_value) except: pass for i in range(len(fields)): # format fields without quotes,dates, or times in them as floats if not any( , 'DATE' in fields , 'TIME' in fields ]): setattr(self, fields , float(values )) else: values = values .replace('', '') setattr(self, fields , values ) # create datetime_obj attribute (drop decimal seconds) dto_string = self.DATE_ACQUIRED + self.SCENE_CENTER_TIME self.DATETIME_OBJ = datetime.strptime(dto_string.split(.) , %Y-%m-%d%H:%M:%S) # only landsat 8 includes sun-earth-distance in MTL file, so calculate it # for the Landsats 4,5,7 using solar module. if not self.SPACECRAFT_ID == LANDSAT_8: # use 0s for lat and lon, sun_earth_distance is not a function of any one location on earth. s = solar(0, 0, self.DATETIME_OBJ, 0) self.EARTH_SUN_DISTANCE = s.get_rad_vector() print(Scene {0} center time is {1}.format(self.LANDSAT_SCENE_ID, self.DATETIME_OBJ))
标签: Landsat OLI
