s5p - Keep Calm

tropomi tools
HARP software can be downloaded from GitHub: download
HARP论坛
NO2 Mapping description, here
MATLAB读取nc文件并转换为.tif格式(here), Matlab：Save as GeoTiff Format (here)
Sentinel-5P数据的栅格化处理 (here)
sp5 harp

编写trans.sh，实现harp和gdal结合的批量nc-nc-tif
conda install -c stcorp harp -y
sudo apt install harp
conda install gdal -y
# harpconvert -a 'keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);bin_spatial(2400,-37,0.01,3300,-83,0.01);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' no2.nc xx.nc

harpconvert -a 'keep(latitude_bounds,longitude_bounds,nitrogendioxide_tropospheric_column);bin_spatial(2400,-37,0.01,3300,-83,0.01);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' no2.nc xx1.nc

全球经纬度的取值范围为：纬度-90至90，经度-180至180
中国的经纬度范围大约为：纬度3.86至53.55，经度73.66至135.05
纬度1度 = 大约111km 
纬度1分 = 大约1.85km 
# harpconvert -a 'keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);bin_spatial(52,3,1,64,73,1);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' no2.nc china.nc
	
坐标换算
# lat_edge_length, lat_edge_offset, lat_edge_step, lon_edge_length, lon_edge_offset, and lon_edge_step.
bin_spatial(181,-90,1,361,-180,1):  latitude starts at -90 with 1 degree steps (ending up after (181 - 1) steps at 90) and longitude starts at -180 with 1 degree steps (ending up after (361 - 1) steps at 180)

1 degree

# harpconvert -a 'keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);bin_spatial(181,-90,1,361,-180,1);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' no2.nc xx1.nc

0.1 degree, 6'=10km

# harpconvert -a 'keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);bin_spatial(1801,-90,0.1,3601,-180,0.1);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' no2.nc f6.nc

0.05 degree, 3'=5km
# harpconvert -a 'keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);bin_spatial(3601,-90,0.05,7201,-180,0.05);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' no2.nc f3.nc

# harpconvert -a 'keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);bin_spatial(3601,-90,0.05,7201,-180,0.05);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' S5P_OFFL_L2__NO2____20200301T053638_20200301T071808_12341_01_010302_20200304T104828.nc tile5f3.nc

# harpconvert -a 'keep(latitude_bounds,longitude_bounds,tropospheric_NO2_column_number_density);bin_spatial(3601,-90,0.05,7201,-180,0.05);squash(time, (latitude_bounds,longitude_bounds));derive(latitude {latitude});derive(longitude {longitude});exclude(latitude_bounds,longitude_bounds,count,weight)' S5P_OFFL_L2__NO2____20200301T071808_20200301T085939_12342_01_010302_20200304T121216.nc tile7f3.nc

# translate to GeoTiff
gdal_translate -a_srs EPSG:4326 -of GTiff NETCDF:"C:\Users\liupei\Desktop\no2\2019\S5P_OFFL_L2__NO2____20190301T031704_20190301T045834_07147_01_010202_20190307T050828.nc.nc":tropospheric_NO2_column_number_density C:/Users/liupei/Desktop/no2/2019/hello.tif
	
# tfw from ipython 学习 gdal https://wenku.baidu.com/view/f419727d27284b73f2425035.html
import gdal
dataset = gdal.Open("f3_S5P_OFFL_L2__NO2____20190301T031704_20190301T045834_07147_01_010202_20190307T050828.nc.tif")

gt = dataset.GetGeoTransform()
gt
gt[1]
gt[2]
gt[4]
gt[5]
gt[0]
gt[3]

dataset.RasterCount
band = dataset.GetRasterBand(1)
dataset.RasterXSize
dataset.RasterYSize

dataset.ReadAsArray(500,500,10,10)
dir(band)
band.XSize
band.YSize
band.DataType
	
import gdalconst
dir(gdalconst)

band.GetNoDataValue()
band.GetMaximum()
band.GetMinimum()
band.ComputeRasterMinMax()

band.GetRasterColorInterpretation()
gdalconst.GCI_PaletteIndex
colormap = band.GetRasterColorTable()
dir(colormap)