Python uses OpenCV library for simple meteorological remote sensing image computing.

Python uses OpenCV library for simple meteorological remote sensing image computing.

The full name of OpenCV is Open Source Computer Vision Library, which is a cross-platform Computer Vision Library. OpenCV is a product launched and developed by Intel. It is issued with BSD license and can be used free of charge in the commercial and research fields. OpenCV can be used to develop real-time image processing, computer vision, and pattern recognition programs. This library can also be accelerated using Intel's IPP.
OpenCV is written in C ++. Its main interface is C ++, but a large number of C language interfaces are retained. This library also has a large number of Python, Java and MATLAB/Ave ave (Version 2.5) interfaces. API interface functions in these languages can be obtained through online documentation. Currently, C #, Ch, and Ruby are also supported.
When compiling OpenCV related to camera input on Windows, some basic classes in DirectShow SDK are required. This SDK can be obtained from the Samples \ Multimedia \ DirectShow \ BaseClasses subdirectory of the pre-compiled Microsoft Platform SDK (or DirectX SDK 8.0 to 9.0c/DirectX Media SDK prior to 6.0.

Next, let's take a look at the application of OpenCV in Python programming. Let's take a look at simple Meteorological Computing and use the opencv library in python to write a script to calculate the image reflectivity of batch processing ~

The core step is remote sensing image spectral radiation calibration, atmospheric correction, and reflectivity calculation.

1. Spectral Radiation Calibration of Remote Sensing Images
Radiation distortion caused by the sensitivity characteristics of a remote sensor is mainly formed by the characteristics of its optical system or photoelectric conversion system. The sensitivity characteristics of the photoelectric conversion system are usually repeated, the calibration is generally carried out through regular ground measurement values.
The following conversion formula is used for remote sensor spectral radiation calibration:

After finding the offset and gain value for each band of the remote sensor, find this table ~

So this function can be scaled:

def computL(gain,Dn,bias):     return (gain*Dn+bias) 

2. Atmospheric correction of Remote Sensing Images
Any atmospheric calibration method that relies on the atmospheric physical model must first perform remote sensor radiation calibration.
This is the formula (Chavez p s, Jr. Image-Based Atmospheric Correction Revisited and Improved Photogrammetric Engineering and Remote Sensing, 1025, 1036)

Among them: Lhazel -- atmospheric spectral radiation value; LI, min -- Minimum spectral radiation value of each remote sensor band; LI, 1% -- radiation value with a reflectivity of 1%.

The formula for the calculation of LI, min, and LI, 1% is omitted. If you are interested, you can check the paper by yourself ~

The parameters required for calculating Lhazel can be obtained from the header file of the remote sensing image, and some are fixed parameters ~ These are all hidden behind the ENVI, but when I write my own scripts, I find that they are still useless.

The code for calculating Lhazel is as follows:

#ESUN ESUNI71=196.9 cos=math.cos(math.radians(90-41.3509605)) # 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 

 

3. Calculate the reflectivity of Remote Sensing Images
Based on the principle and process of solar radiation and atmospheric transmission, the mathematical model of the TM/ETM + Data land reflectivity inversion can be comprehensively expressed:

Among them: p -- relative reflectivity of the ground; D -- astronomical unit distance of the day and ground; LsatI -- spectral radiation value of the sensor, that is, the radiation energy of the top layer of the atmosphere; LhazeI -- atmospheric radiation value; ESUNl: average spectral radiation of the top layer of the atmosphere, that is, the solar radiation of the top layer of the atmosphere; SZ: the top corner of the sun.

The formula for calculating the two parameters is as follows:
Day-to-earth astronomical unit distance D = 1-0.01674 cos (0.9856 × (JD-4) × π/180 );
(JD is the Julian Day of Remote Sensing Imaging. The formula is as follows:

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 

I, J, and K are respectively years, months, and days.

With this, we can finally calculate the reflectivity directly. The rough code is as follows, because it is written to play, and there is no way to deal with it:
However, when using remote sensing images for computing and output, you need to use a uint16 array (uint8 is not long enough ..)
Some parameters involve floating-point number calculation. If there are extremely high requirements on the processing results, it is best to use a dedicated scientific computing Library (such as me, it doesn't matter)

Import cv2 import numpy as np import math img1 = cv2.imread ('f: \ L71121040_04020030220_B10.TIF ') # image format conversion img10 = cv2.cvtColor (img1, cv2.COLOR _ BGR2GRAY) # computing 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 # Set ESUNI value ESUNI71 = 196.9 # Calculate the daily distance from D = 1-0.01674 * math. cos (0.9856 * (JD-4) * math. pi/180) # Calculate the top corner of the Sun cos = math. cos (math. radians (90-41.3509605) inter = (math. pi * D)/(ESUNI71 * cos) # atmospheric calibration parameter settings Lmini =-6.2 Lmax = 293.7 Qcal = 1 Qmax = 255 LIMIN = Lmini + (Qcal * (Lmax-Lmini)/Qmax) LI = (0.01 * ESUNI71 * cos)/(math. pi * 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 + bias) if _ name _ = '_ main _': print 'd = ', D print 'coszs =', cos print 'lhazel = ', lhazel # Calculate image reflectivity result = np. zeros (img. shape, dtype = 'int16') for I in range (0, img. shape (1): for j in range (0, img. shape (0): Lsat = computL (1.18070871, img10 [I, j],-7.38070852) result [I, j] = inter * (Lsat-Lhazel) * 1000 # Save the image cv2.imwrite ("F: \ result. tif ", result) cv2.namedWindow (" Image ") cv2.imshow (" Image ", result) cv2.waitKey (0)

Articles you may be interested in:

• Python and opencv for Face Detection and Tracking
• How to Use OpenCV to rotate images in Python
• Tutorial on installing Python and OpenCV on Raspberry Pi 2 or Raspberry Pi B +
• Using Python and OpenCV libraries to convert URLs to OpenCV formats
• An example of using OpenCV for Face Detection in python