Using the OpenCV library for simple meteorological remote sensing image calculation in Python

The full name of OpenCV is the open Source Computer Vision Library, which is a cross-platform computer vision database. OPENCV is initiated and developed by Intel Corporation and is licensed for distribution in the BSD license and is available free of charge in the commercial and research fields. The OPENCV can be used to develop real-time image processing, computer vision, and pattern recognition programs. The library can also be accelerated with Intel Corporation's IPP.
OpenCV is written in the C + + language, and its main interface is the C + + language, but still retains a large number of C-language interfaces. The library also has a large number of interfaces for Python, Java and Matlab/octave (version 2.5). API interface functions for these languages are available through online documentation. Support for C #, Ch,ruby is now also available.
Some of the base classes in the SDK need to be DirectShow when compiling OPENCV on Windows with the relevant parts of the camera input. The SDK can be from a subdirectory of the precompiled Microsoft Platform SDK (or DirectX SDK 8.0 to 9.0c/directx Media SDK prior to 6.0) Samples\multimedia\dir Ectshow\baseclasses obtained.

Let's take a look at the application of OpenCV in Python programming, let's deal with simple meteorological calculations, write a script using the OpenCV library in Python to calculate the reflectivity of a batch image, and try it out.

The core step is the remote sensing image spectral radiometric calibration → atmospheric correction → calculate reflectivity This is three steps.

1. Spectral radiometric calibration of remote sensing imagery
The radiation distortion caused by the sensitivity characteristic of the remote sensor is mainly formed by the characteristic of its optical system or photoelectric conversion system, and the sensitivity characteristic of photoelectric conversion system is often repeated, and its correction is usually carried out by the regular ground measurement value.
The following conversion equations are used for the spectral radiometric calibration of remote sensors:

The shift and gain value of each band of the remote sensor from the paper find out, find such a table ~

So this function can be fixed:

def computl (Gain,dn,bias):     

2. Atmospheric correction of remote sensing imagery
Any method of atmospheric calibration that relies on atmospheric physical models requires radiometric calibration of the remote sensor first.
The formula is this slightly (Chavez P s,jr. image-based atmospheric Correction revisited and improved Photogrammetric Engineering and Remote Se Nsing, 1996,62,1025-1036)

Wherein: lhazel--atmospheric spectral radiation value; The minimum spectral radiation value of each band of li,min--remote sensor; li,1%--the blackbody radiation value with a reflectance of 1%.

On the li,min and li,1% the calculation formula is omitted Ah, interested students can check their own papers ~

The parameters needed to calculate the Lhazel can be obtained from the remote sensing image header file, and some are fixed parameters ~ These are hidden behind the envi, but when they write the script to find out they are still a waste of kung fu.

The code to calculate the Lhazel is as follows:

3. Calculate the reflectance of remote sensing images
According to the principle and process of solar radiation and atmospheric transmission, the mathematical model of ground reflectance inversion of tm/etm+ data can be expressed as:

Among them: ρ--ground relative reflectance, d--astronomical unit distance; lsati--sensor spectral radiation value, that is, the radiant energy of the atmospheric top layer; lhazei--atmospheric radiation value; The solar average spectral radiation of the esunl--atmosphere, i.e. the solar irradiance of the atmospheric top layer; sz--Sun Zenith Angle.

Here is a calculation formula for two of these parameters:
The Astronomical unit distance d=1-0.01674 cos (0.9856x (JD-4) xπ/180);
JD is a Julian day for remote sensing imagery, calculated as:

I, J, K respectively for the year, month, day

With these, the final can be directly calculated the reflectivity, rough code as follows, because it is written to play, nor how to deal with:
However, it is important to note that the remote sensing image is computed and output, and needs to be stored using an array of type uint16 (uint8 length is not enough.) ）
Some parameters are related to floating-point calculation, if the processing results are very high, it is best to use a specialized scientific database (like my slag school do not mind these)

Import CV2 Import NumPy as NP import math Img1=cv2.imread (' F:\L71121040_04020030220_B10. TIF ') #图像格式转换 Img10=cv2.cvtcolor (Img1,cv2. Color_bgr2gray) #计算JD i=2003 j=2 k=20 jd=k-32075+1461* (i+4800+ (J-14)/12)/4+367* (j-2-(J-14)/12*12)/12-3* ((I+4900+ ( J-14) (/12)/100)/4 #设置ESUNI值 esuni71=196.9 #计算日地距离D D=1-0.01674*math.cos ((0.9856* (JD-4) *math.pi/180) #计算太阳天顶角 cos= Math.Cos (Math.radians (90-41.3509605)) inter= (math.pi*d*d)/(Esuni71*cos*cos) #大气校正参数设置 lmini=-6.2 Lmax=293.7 Qcal=1 qmax=255 limin=lmini+ (qcal* (Lmax-lmini)/qmax) li= (0.01*esuni71*cos*cos)/(math.pi*d*d) lhazel=limin-li def copy (IMG, New1): new1= np.zeros (img.shape,dtype= ' uint16 ') new1[:,:] = img[:,:] def computl (Gain,dn,bias): Return (Gain*dn+bia s) If __name__ = = ' __main__ ': print ' d= ', D print ' coszs= ', cos print ' lhazel= ', Lhazel #计算图像反射率 Result=np.zeros ( Img.shape,dtype= ' uint16 ') for I in Range (0,img.shape (1)): for J in Range (0,img.shape (0)): Lsat=computl (1.18070 871,img10[i,j],-7.38070852) result[i,j]=inter* (Lsat-lhazel) *1000 #保存图像 cv2.imwrite ("f:\\result.tif", result) Cv2.namedwindow ("Image") cv2.imshow  ("Image", result) cv2.waitkey (0)

