![](https://img-blog.csdnimg.cn/20201014180756923.png?x-oss-process=image/resize,m_fixed,h_64,w_64)
遥感
文章平均质量分 55
江北20190411
。
展开
-
Google Earth Engine (GEE) 提取某一位置时间序列值
前言做时间序列相关算法的同学,经常需要下载年际的数据来进行试验。而且算法研究的比较重要的步骤就是利用实测站点数据进行验证。本文讲述如何使用GEE下载指定卫星遥感/再分析数据集指定位置(实测站点位置)的时间序列值。提取指定点年序列值-imageCollection整个过程可以在imageCollection上进行操作,并且中间处理过程的返回对象都是image,这种情景非常简单直接。配合上篇博客,也可以进行QC过滤。首先定义点的经纬度坐标。然后过滤数据集合,这里以MODIS的NDVI产品MOD13A2为原创 2022-03-20 17:55:25 · 6994 阅读 · 5 评论 -
Google Earth Engine (GEE) 实现对MODIS产品批量质量控制
前言最近处在学位论文初稿完成和申博之间的空档期,所以想学点新东西。GEE火了很久了,一直没用真正使用过,还是坚持下载数据到本地,然后用Python处理。主要是怕一旦接触GEE,就会花费很多心思在上面。再者,我本身是做算法的,不仅仅是数据分析,需要使用本地数据来应用算法。现在觉得下载到本地然后处理数据有很多不确定性,不如直接在GEE上完成预处理然后导出,可以省多很多精力给算法研究本身。所以,今后主要关注于使用GEE进行数据预处理和导出的相关实现。言归正传,https://spatialthoughts.c原创 2022-03-13 22:04:48 · 2951 阅读 · 4 评论 -
Python GDAL工具使用及使用VRT格式数据处理
前言VRT及虚拟栅格文件。使用方法先找到gadl工具exe文件的路径D:\Anaconda\envs\geopandas\Library\bin对数据列表构建vrt:首先在工具路径下打开cmd,或者打开cmd输入完整工具路径构建虚拟栅格gdalbuildvrt -input_file_list E:\GeoLearn\dem_vrt_test\data_list.txt E:\GeoLearn\dem_vrt_test\data_vrt.vrt对虚拟栅格进行重投影, 并指定输出分辨率原创 2021-09-28 11:26:46 · 2637 阅读 · 0 评论 -
Python计算灰度共生矩阵和常用的纹理因子
前言如何使用skimage计算GLCM和纹理因子可以参考利用python的skimage计算灰度共生矩阵。但是Skimage库只提供了常用的8中纹理因子(均值、方差、同质性、对比度、差异性、熵、角二阶矩、相关性、)中的5种,缺少均值、方差和熵的计算。这里讲解如何修改源码添加这三个因子的计算。添加后调用方式和其他因子的计算相同。修改/添加部分找到 skimage库greycoprops源码。Spyder可以直接右键该函数,转到定义处。需要添加/修改的部分,保持了与源码相同的编码风格,全是numpy向量原创 2021-09-28 11:19:06 · 2014 阅读 · 7 评论 -
累计分布函数匹配 (CDF matching) Python实现
import numpy as npimport scipy.io as ioimport matplotlib.pyplot as pltfrom scipy.interpolate import InterpolatedUnivariateSplinefrom scipy.optimize import curve_fitnp.random.seed(0)# 原始数据x = np.arange(0, 37)data_or = np.sin(x) + np.random.normal(s原创 2021-09-28 10:39:28 · 1588 阅读 · 0 评论 -
ASTER产品数据 (AST_08) 几何校正
前言国内有关ASTER的博客等资料真的很少。因个人需求,用过这个数据的不少产品,分享一下校正的方法。数据可从[EarthData](https://search.earthdata.nasa.gov/)进行下载,需要注册账号。原创 2021-05-17 16:37:30 · 1570 阅读 · 5 评论 -
Python scipy数组线性插值
实验数据是一个35,35的数组,采用线性插值得到375,375的数组。构造interp2d对象时,要求输入原始数据的行列坐标和数值,并指定插值方式。然后向构造好的对象中输入想要得到的数组坐标值。import matplotlib.pyplot as pltimport numpy as npfrom scipy.interpolate import interp2dresidual = np.load(r'F:\Py_project\Spatial_RF\residual.npy').res.原创 2021-05-17 15:35:17 · 841 阅读 · 0 评论 -
FY-4A 静止卫星圆盘数据几何校正
前言最近好多网友有FY-4的几何校正需求,就来更新一下。数据准备L1级的FY-4数据是HDF格式的全圆盘数据,数据集里面没有经纬度信息。每个位置对应的经纬度数据需要单独下载,如下图。这里我下载了4 km的经纬度查找表。文件下载解压后,是一个.raw格式的二进制文件。下面用python解析这个文件。import numpy as npfile = r'F:\Py_project\geoFY-4A\FullMask_Grid_4000.raw'raw_image = np.fromfile原创 2021-04-20 15:21:18 · 4601 阅读 · 9 评论 -
Python+Cartopy绘制中国地图
前言之前看了很多cartopy的教程,刚好手上有一份很久前下载的ERA5水汽数据,就练习一下吧。maskout白化的代码和colorbar分割来自气象家园的大佬。代码及数据网盘链接提取码:vtmk代码import cartopy.crs as ccrsimport matplotlib.pyplot as plt import numpy as npimport netCDF4 as ncimport cartopy.feature as cfeatureimport cartopy.io原创 2021-04-18 11:46:13 · 5402 阅读 · 6 评论 -
Python绘制分类图
前言遥感影像分类图一般为特定数值对应一类地物,用Python绘制时,主要在颜色的映射和对应的图例生成。plt.matplotlib.colors.ListedColormap支持自定义颜色。matplotlib.patches mpatches对象可以生成一个矩形对象,控制其颜色和地物类型的颜色对应就可以生成地物分类的图例了。具体用法可以自行Google和百度。下面给出一个模拟地物分类数据的可视化例子。代码import numpy as npimport matplotlib.pyplot as p原创 2021-04-17 18:03:35 · 1459 阅读 · 4 评论 -
Python+Basemap 经纬度标签的旋转与间距调整
前言这是个非常个人的需求,原图如下,可以看到经度的标注是水平的,多子图并列的时候占空间很大,也不够美观。所以目的是将它垂直显示。垂直显示经度标签basemap绘制投影影像参看前面的博客。下面两句是分别是绘制经纬线,查看它俩的方法和属性,发现并没有控制标签角度的部分。但是basemap是基于matplotlib库的,所以可以通过操作matplotlib来实现目的。drawparallels = m.drawparallels(parallels,labels=[True,False,False,原创 2021-03-15 21:11:55 · 1708 阅读 · 0 评论 -
Py6s批量大气校正Landsat8
前言使用py6s进行Landsat8等遥感影像的大气校正。代码来自基于6s模型的遥感影像大气校正方法,对应的代码下载地址Github。由于原始代码是在cmd中通过命令行传参运行的,包括两个输入参数,分别是 Input_dir=输入路径 Output_dir=输出路径。这样每次只能处理一个时相的数据。本文目的在于提供一个Windows环境下通过不修改源码的形式来实现批量大气校正多时相影像的思路。处理单时相将Github上整个代码工程下载到本地。打开Anaconda Prompt,输入python "原创 2021-03-14 09:20:20 · 3389 阅读 · 13 评论 -
Python+Cartopy绘制已投影影像
前言数据是一副投影为UTM 50°N的地温影像, 使用ax.contourf绘制已投影栅格时需要指定参考坐标系ccrs 和每一个数据点的投影坐标系坐标。(2) 数据读取部分的代码省略,是通过gdal读取的,readTifAsArray函数前面的博客中可以找到。代码指定投影并绘制import cartopy.crs as ccrsfrom cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter, LatitudeLocator原创 2021-02-24 09:46:34 · 1651 阅读 · 2 评论 -
SURFRAD (Surface Radiation Budget) Network 数据下载
下载地址Global Monitoring Laboratory点进任意一个站点FTP,很清楚的看到数据的下载列表和READMEREADME对数据质量的说明:文件第138行给出了数据的内容和格式的描述需要注意的是每一种实测数据成列分布,质量控制代码紧跟在实测数据之后,如下所示。计算地表温度所需的F↓和F↑分别对应数据中以下内容。...原创 2020-10-12 11:03:24 · 1665 阅读 · 2 评论 -
Python+Basemap绘制已投影的影像
前言数据一副具有albers投影的地温影像,而Basemap的基本参数中要求输入影像的左下右上的经纬度坐标,所以关键在于如何将投影坐标转化为大地经纬度坐标。在坐标转换过程中,用到了pyproj库,首先要定义转换前后的坐标类型,定义WGS84坐标系可用pyproj.CRS.from_epsg(4326)进行定义,但是由于我的pyproj的数据库路径有些问题,所以首先采用gdal读取一幅WGS84坐标的影像的wkt投影信息,然后pyproj.CRS.from_wkt( )进行定义要转换成的坐标系。关于影像原原创 2020-09-27 09:54:24 · 1933 阅读 · 0 评论 -
遥感辐射亮度单位转换
关于波数与波长之间的转换参考如下:these two units do not have a linear scaling relationship. rather, this depends on wavelength. wavenumber ν [cm-1] relates to wavelength λ [nm] as v = 10^7/λ. when integrating the same radiances (in their appropriate units) over the same原创 2020-09-18 09:54:07 · 1518 阅读 · 2 评论 -
Python+GDAL几何校正任意自带经纬度数据的遥感影像
# -*- coding: utf-8 -*-import gdal, osrimport os def geoMERSI2(file, geoFile, afterGeoFile): ''' Parameters ---------- file : 文件绝对路径 需要校正的文件. geoFile : 文件绝对路径 地理坐标存在的文件. afterGeoFile : tif文件绝对路径 经过校正原创 2020-09-12 15:41:31 · 5209 阅读 · 25 评论 -
Python实现多探元并扫遥感影像的直方图匹配
前言直方图匹配的函数参考了https://www.jianshu.com/p/3f6abf3eeba2写的很好。import matplotlib.pyplot as pltimport numpy as npfrom PIL import Imageimport matplotlibimport gdal, osrfrom osgeo import gdal_arrayimport pandas as pd# 将读取文件的灰度矩阵,转化为直方图,这里的直方图定义为python的dict类原创 2020-09-10 12:33:25 · 1041 阅读 · 2 评论 -
ASTER光谱库批量转换为传感器通道比辐射率
前言最近做的一个实现。将ASTER光谱库的所有地物连续的波谱响应的制定地物大类向指定卫星的通道比辐射率进行转换。最后形式为pd.Dataframe,列索引为地物名字,值为通道比辐射率。import osimport numpy as npfrom shutil import *import matplotlib.pyplot as plt from scipy import integrateimport pandas as pd#获取指定路径下所有文件的绝对路径def listDir(原创 2020-09-10 12:24:15 · 675 阅读 · 2 评论 -
Python读取nc文件转tif
import netCDF4 as ncdata = r'Y:\folder\H08_20200601_0000_1D_ROC010_FLDK.02401_02401.nc'nc_data_obj = nc.Dataset(data)nc_data_objOut[4]: <class 'netCDF4._netCDF4.Dataset'>root group (NETCDF4 data model, file format HDF5): title: Himawari-原创 2020-09-10 12:13:53 · 9370 阅读 · 10 评论 -
TIGR2000大气廓线数据库下载地址
The Thermodynamic Initial Guess Retrieval (TIGR) database TIGR2002,constructed by the Laboratoire de Meteorologie Dynamique and representing a worldwide set of atmospheric situations (2311 radisoundings) from polar to tropical atmospheres with varying wat.原创 2020-08-04 16:01:30 · 2179 阅读 · 7 评论 -
ENVI图像自动配准工具(转载)
ENVI 5.0增加了Image Registration Workflow,它是一个全新的影像配准工作流,具有自动、准确、快速的特点。它将之前版本中繁杂的参数设置步骤集成到统一的面板中,并且增加了生成种子点的影像匹配参数设置项、Harris角点算子、匹配粗差踢除算子。在少量或者无需人工干预的情况下该影像配准工作流能快速而准确的实现影像间的自动配准。功能详细操作过程 启动ENVI5,在ToolBox中选择GeometricCorrection->Regi转载 2020-06-02 18:09:47 · 1753 阅读 · 0 评论 -
MODIS产品质量控制文件使用方法
官方关于产品质量控制的说明(机翻)质量指标 在生产过程中生成的CoreMetadata.0全局属性QA 中的元数据对象以及质量控制(QC)SDS中给出,或者在数据产品的产品后科学和质量检查中给出。CoreMetadata.0中的 QA元数据对象全局属性是AutomaticQualityFlag和ScienceQualityFlag及其相应的说明。根据运行LST算法期间遇到的数据条件,根据规则设置AutomaticQualityFlag。此质量检查标志的设置是完全自动化的。设置它的规则是自由的;几乎所有数原创 2020-05-18 00:51:13 · 6912 阅读 · 13 评论 -
IDL自定义滑动窗口函数
许多地表参数反演算法中是以整景影像中的每个小窗口为单位进行计算。这时候就需要对整景影像划分为数个固定窗口大小影像。以下函数实现输入想要设定窗口的大小,影像列行数,返回一个分块影像的左上角和右下角坐标的查找表数组。有这样一个数组,无论是影像计算还是裁剪深度学习样本都可以适用。原创 2020-05-11 23:04:10 · 1346 阅读 · 0 评论 -
ASTER连续光谱向传感器通道光谱转换
数据说明 ASTER光谱库包含三个子库,分别是 Johns Hopkins University (JHU) 光谱库、Jet Propulsion Laboratory (JPL) 光谱库和 United States Geological Survey (USGS) 光谱库。ASTER 光谱库包含多达 2300 余种典型地物的光谱,包括矿物质、岩石、土壤、植被、水体、冰雪等各类自然地物以及人造地物等类别,其波谱范围从 0.4 μm 到 15.4 μm。问题描述 但由原创 2020-05-09 16:17:46 · 1555 阅读 · 13 评论 -
Python GDAL批处理MODIS产品
批量几何校正MODIS产品from osgeo import gdal, osrimport osimport time'''批量读取文件夹中hdf格式的modis产品数据并进行几何校正(WGS84)'''def readHdfWithGeo(hdfFloder, saveFloder): #获取输入文件夹中所有的文件名 hdfNameList = os.listdi...原创 2020-05-07 22:10:56 · 2735 阅读 · 2 评论 -
Python GDAL学习笔记(二)
目录一、自定义空间坐标系使用EPSG或EPSGA编号进行初始化使用WKT字符串进行初始化使用PROJ字符串进行初始化附录二、遥感影像几何校正利用卫星自带的地理定位文件进行几何校正查看数据自带经纬度查找表位置查找需要进行校正的数据进行地理校正三、空间参考系统转换一、自定义空间坐标系使用EPSG或EPSGA编号进行初始化EPSG:欧洲石油勘探组织(European Petroleum Surve...原创 2020-05-05 17:07:50 · 6930 阅读 · 7 评论 -
Python GDAL学习笔记(一)
目录导入GDAL模块并查看版本/位置打开影像查看行列号和波段数查看影像文件描述和其他属性打开数据集波段并获取统计信息利用numpy将波段数据写入数组将波段数据写入数组的两种方法对比导入GDAL模块并查看版本/位置from osgeo import gdalprint("GDAL's version is:" + gdal.__version__)print(gdal)GDAL's ve...原创 2020-05-04 17:05:10 · 7645 阅读 · 3 评论 -
矩匹配算法影像去噪
2020-5-2-9:19常规的矩匹配算法能起到去除条带噪声的效果,但去噪后影像灰度信息损失过大。IDL关键代码: dim = size(data, /dimensions) row = dim[1] col = dim[0] re_mean = mean(mean(data, dimension=1)) re_stev = mean(stddev(data, dimension=1)) or_mean = mean(data, dimension=1) or_st原创 2020-05-03 21:59:05 · 2627 阅读 · 1 评论 -
像元二分模型计算植被覆盖度
前言关于改进的像元二分模型第一次出现是在李苗苗老师的博士论文中,但她后来发表的一篇期刊文章把这一部分单独摘了出来并做了应用案例,期刊链接在下方。密云水库上游植被覆盖度的遥感估算论文要点计算公式为:Fc =(NDVI -NDVIsoil) (NDVIveg -NDVIsoil)问题的关键在于NDVIveg 与NDVIsoil如何取值。 其中 ,NDVIsoil 为完全是裸土或无植...原创 2020-05-01 00:43:52 · 11396 阅读 · 0 评论 -
HWSD全球土壤数据下载处理
下载地址:http://www.fao.org/soils-portal/soil-survey/soil-maps-and-databases/harmonized-world-soil-database-v12/en/原创 2020-04-30 10:33:57 · 27284 阅读 · 40 评论 -
6S模型逐像元大气校正MERSI2
前言FY3D-MERSI2 6S大气校正IDL调用6S模型之前完成了6S模型源码中添加MERSI2的光谱响应函数和实现了IDL调用6S模型,然后进行了应用同一参数进行大气校正。这次实现利用IDL调用6S模型建立大气校正系数查找表然后对MERSI2进行逐像元大气校正。由于MERSI2自带的角度信息是逐像元的数据,而我建立的查找表步长比较大,因此在进行校正前,先对角度信息和MODIS AOD...原创 2020-04-28 18:15:34 · 1480 阅读 · 0 评论 -
FY3D MERSI2 6S大气校正
终于在6S源码里添加了MERSI2 650、865nm波段的响应函数。网上的博客已经写的很详细了,需要注意的是更改源码的时候最好在编译器里进行更改。我一开始是在NotepadC++里更改的,然后到Visual Fortran 6.6里进行编译的时候,各种变量,格式报错。之后问了同门,他给我推荐了Code::Blocks,最后终于顺利编译。注意最后编译的时候,6S源码的路径不要太深,也要避免中文...原创 2020-04-25 19:47:27 · 1890 阅读 · 1 评论 -
传感器通道波长单位换算
基础概念波数(k):为波长λ的倒数,即1cm中所含波的个数。原子、分子和原子核的光谱学中的频率单位.符号为σ或v.等于真实频率除以光速,即波长(λ)的倒数,或在光的传播方向上每单位长度内的光波数.其常用单位为cm-1,SI制单位为m-1.波长(λ):是指波在一个振动周期内传播的距离。也就是沿着波的传播方向,相邻两个振动位相相差2π的点之间的距离。换算关系波数(cm-1) ...原创 2020-04-23 20:05:06 · 1501 阅读 · 1 评论 -
IDL矢量/掩膜裁剪影像(ENVI API编程)
前言当我们想要用一个矢量文件去裁剪影像时,可以使用ENVIVectorMaskRaster函数。maskedRaster = ENVIVectorMaskRaster(raster, shp_file)**图中标记红圈处是我们裁剪出的影像,非常小的一块。但是我却发现用ENVI-IDL裁剪出的影像仍然保留着原来影像的行列数,并不像ENVI中Subset Data from ROIS工具一...原创 2020-04-22 21:13:02 · 5411 阅读 · 6 评论 -
IDL实现FY3D MERSI2 250M/1KM数据预处理
结合我自己的需求只处理了Band3.4.15.18.24.25,后续获取光谱响应函数后会加上6s大气校正部分。由于构建250M的GLT文件时间比较长,所以没有添加构建GLT部分,构建GLT的方法可以参考ENVI-IDL技术殿堂官方博客。代码末删除临时文件的方法比较暴力,带有[envitempfile]的文件都会删除,结合实际情况,谨慎使用!原创 2020-04-19 22:09:00 · 3492 阅读 · 6 评论 -
ENVI下预处理GF2-PMS1数据
数据数据是一景2018年11月27日GF2-PMS1-L1A级影像。处理步骤正射校正辐射定标、大气校正FLAASH大气校正设置如下:处理结果大气校正前后植物光谱曲线对比:融合前后影像对比:...转载 2020-04-17 16:35:39 · 2275 阅读 · 0 评论 -
哨兵2号真彩色影像合成
所需软件SNAP下载地址 Sen2Cor下载地址数据下载欧空局下载地址处理过程一、L1C级产品大气校正(1)在Sen2Cor的安装目录的地址栏输入cmd打开命令窗口(2)在cmd窗口中输入L2A_Process.bat --resolution 10 E:\snap\S2AL1C\S2A_MSIL1C_20200130T025941_N0208_R032_T50SNF_2020...原创 2020-04-17 00:49:47 · 9457 阅读 · 3 评论 -
FY3D MERSI2定标问题
问题描述首先根据手册定标Band3得到ref。然后查看影像值并圈定一个roi进行统计得到如下:此时发现影像值并不介于0-1,但手册明确说明定标后为反射率ref。于是猜测是否还要除以100来转换到0-1之间的反射率。验证验证的方法是选取与MERSI2过境时间相近的Modis影像,然后将MERSI2的辐亮度与MODIS的辐亮度进行对比。首先根据如下公式,将MERSI2 band3定标后...原创 2020-04-18 19:26:43 · 2030 阅读 · 11 评论 -
遥感物理笔记
辐射量度概念辐射能量(W):电磁辐射的能量,单位:J;辐射通量Φ;单位时间内通过某一面积的辐射能量,Φ=dW/dt,单位是w;辐射通量是波长的函数,总辐射通量应该是各波段辐射通量之和或辐射通量的积分值。辐射通量密度(E):单位时间内通过单位面积的辐射能量,E=dΦ/dS,单位:w/m2,S为面积。辐射照度(I):被辐射的物体表面单位面积上的辐射通量,I=dΦ/dS,单位是w/m2, s为面...原创 2020-04-04 23:00:23 · 4215 阅读 · 0 评论