江北20190411-CSDN博客

原创 Google Earth Engine (GEE) 提取某一位置时间序列值

前言做时间序列相关算法的同学，经常需要下载年际的数据来进行试验。而且算法研究的比较重要的步骤就是利用实测站点数据进行验证。本文讲述如何使用GEE下载指定卫星遥感/再分析数据集指定位置（实测站点位置）的时间序列值。提取指定点年序列值-imageCollection整个过程可以在imageCollection上进行操作，并且中间处理过程的返回对象都是image，这种情景非常简单直接。配合上篇博客，也可以进行QC过滤。首先定义点的经纬度坐标。然后过滤数据集合，这里以MODIS的NDVI产品MOD13A2为

2022-03-20 17:55:25 7091 5

原创 Google Earth Engine (GEE) 实现对MODIS产品批量质量控制

前言最近处在学位论文初稿完成和申博之间的空档期，所以想学点新东西。GEE火了很久了，一直没用真正使用过，还是坚持下载数据到本地，然后用Python处理。主要是怕一旦接触GEE，就会花费很多心思在上面。再者，我本身是做算法的，不仅仅是数据分析，需要使用本地数据来应用算法。现在觉得下载到本地然后处理数据有很多不确定性，不如直接在GEE上完成预处理然后导出，可以省多很多精力给算法研究本身。所以，今后主要关注于使用GEE进行数据预处理和导出的相关实现。言归正传，https://spatialthoughts.c

2022-03-13 22:04:48 3029 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 2741

原创 Python计算灰度共生矩阵和常用的纹理因子

前言如何使用skimage计算GLCM和纹理因子可以参考利用python的skimage计算灰度共生矩阵。但是Skimage库只提供了常用的8中纹理因子（均值、方差、同质性、对比度、差异性、熵、角二阶矩、相关性、）中的5种，缺少均值、方差和熵的计算。这里讲解如何修改源码添加这三个因子的计算。添加后调用方式和其他因子的计算相同。修改/添加部分找到 skimage库greycoprops源码。Spyder可以直接右键该函数，转到定义处。需要添加/修改的部分，保持了与源码相同的编码风格，全是numpy向量

2021-09-28 11:19:06 2038 7

原创 Python绘制双变量三维曲面图

import numpy as npimport matplotlib.pyplot as plt# 定义双变量函数def func(X): x1, x2 = X return (x1-3.14)**2 + (x2-2.72)**2 + np.sin(4*x1+1.41) - np.cos(4*x2-1.73)# 生成双变量网格点x = np.arange(10).reshape(-1, 1)x1, x2 = np.meshgrid(x, x)# 求解y = np.as

2021-09-28 11:06:01 1372

原创经验正交函数 (EOF) / 主成分 (PCA) 分解及Python实现

前言EOF 经验正交函数分解原理及用法：https://wenku.baidu.com/view/93bc8a9b680203d8ce2f2468.html方法实现过程EOFPCAPython实现import numpy as npfrom scipy import linalg# 生成数据并打乱np.random.seed(9)data = np.random.randint(low=1, high=50, size=600).reshape(20, 30) * 10#

2021-09-28 11:00:04 3295

原创差分进化算法 (DE) Python实现及源码解读

设函数，有三个未知解def obj_func(p): x1, x2, x3 = p return x1 ** 2 + x2 ** 2 + x3 ** 2初始解形状self.X = None # shape = (size_pop, n_dim)生成初始解（DE继承类GeneticAlgorithmBase）lb, ub为x的上下限，新解为约束条件上下界内均匀分布的向量class DE(GeneticAlgorithmBase): def __init__(self, f

2021-09-28 10:51:33 2247 2

原创累计分布函数匹配 (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 1645

原创 ASTER产品数据 (AST_08) 几何校正

前言国内有关ASTER的博客等资料真的很少。因个人需求，用过这个数据的不少产品，分享一下校正的方法。数据可从[EarthData](https://search.earthdata.nasa.gov/)进行下载，需要注册账号。

2021-05-17 16:37:30 1637 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 844

原创 GWR4 软件输入数据制作

前言GWR软件是实现地理加权回归建模的专业软件。GWR软件下载地址Click here。软件下载后自带了User guidance，也可以参照这篇博客。软件要求输入的数据格式为.csv。数据中要包括投影坐标，因变量，解释变量。数据制作我这里以地表温度(LST)做因变量，NDVI做解释变量为例，使用Arcgis制作.csv输入文件。LST和NDVI的原始数据是两幅投影坐标的栅格影像。首先在Arcgis中打开影像，然后使用栅格转点工具将栅格转换为点状shp文件。得到每个像素点的shp文件。

2021-04-26 10:53:05 2344 8

原创 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 4763 10

原创 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 5463 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 1472 4

原创 Python+Basemap 经纬度标签的旋转与间距调整

前言这是个非常个人的需求，原图如下，可以看到经度的标注是水平的，多子图并列的时候占空间很大，也不够美观。所以目的是将它垂直显示。垂直显示经度标签basemap绘制投影影像参看前面的博客。下面两句是分别是绘制经纬线，查看它俩的方法和属性，发现并没有控制标签角度的部分。但是basemap是基于matplotlib库的，所以可以通过操作matplotlib来实现目的。drawparallels = m.drawparallels(parallels,labels=[True,False,False,

2021-03-15 21:11:55 1747

原创 Py6s批量大气校正Landsat8

前言使用py6s进行Landsat8等遥感影像的大气校正。代码来自基于6s模型的遥感影像大气校正方法，对应的代码下载地址Github。由于原始代码是在cmd中通过命令行传参运行的，包括两个输入参数，分别是 Input_dir=输入路径 Output_dir=输出路径。这样每次只能处理一个时相的数据。本文目的在于提供一个Windows环境下通过不修改源码的形式来实现批量大气校正多时相影像的思路。处理单时相将Github上整个代码工程下载到本地。打开Anaconda Prompt，输入python "

2021-03-14 09:20:20 3492 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 1671 2

原创 SURFRAD (Surface Radiation Budget) Network 数据下载

下载地址Global Monitoring Laboratory点进任意一个站点FTP，很清楚的看到数据的下载列表和READMEREADME对数据质量的说明：文件第138行给出了数据的内容和格式的描述需要注意的是每一种实测数据成列分布，质量控制代码紧跟在实测数据之后，如下所示。计算地表温度所需的F↓和F↑分别对应数据中以下内容。...

2020-10-12 11:03:24 1706 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 1948

原创遥感辐射亮度单位转换

关于波数与波长之间的转换参考如下：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 1545 2

原创解决：object of type ＜class ‘numpy.float64‘＞ cannot be safely interpreted as an intege

回退一下numpy的版本：conda install numpy==1.17https://github.com/tensorflow/models/issues/8576

2020-09-14 14:58:10 897

原创 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 5336 25

原创解决：TypeError: handle_get_file_code() got an unexpected keyword argument ‘save_all‘

问题出现在更新spyder到最新版本conda install spyder=4.1.4spyder更新完成后，之前的代码都不能正常运行报错：This command failed to be executed because an error occurred while trying to get the file code from Spyder's editor. The error was:An exception has occurred, use %tb to see the ful

2020-09-10 16:24:43 8744 6

原创 numpy.where()多条件判断

np.where中多条件使用 & 或者 |原因：& 和 | 可以判断np.array()格式的true与false数组，而and，or 只能判断单个值true或者false如下所示。

2020-09-10 12:42:38 3473

原创 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 1057 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 690 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 9450 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 2228 7

原创 Origin绘制双x,y轴图像

1.前言目的是绘制双x，双y轴的图像，并且两组数据的个数不一致。2.准备模拟数据首先在excel中模拟了两组数据，P列为1-5随机个数的数值，用于绘制直方图。Q，S列的关系为S=sin(Q),用于绘制折线图。3.绘图将P列数据粘贴至Origin，选中整列右键将类型设置为忽略。点击绘图，选择绘制直方图。双击图窗中图像的上坐标轴，将下图的对号取消，即不显示上轴，同样方法取消右轴的显示。接下来在图窗非图像的空白图右击，选择添加如下轴。将另外两列数据添加至Origin，并修改类型分别为

2020-07-24 16:09:03 17963

原创解决Origin导出图片失真问题

导出图片时，如果是png或者tif等格式，色深选择24位彩色。

2020-07-24 15:33:39 12636 4

原创图像客观质量评价指标IDL实现

;用于计算原始影像与经过处理后影像的平均绝对误差 MAEfunction calMAE, image_or, image_deal nDim = size(image_or, /n_dimension) nDim2 = size(image_deal, /n_dimension) if nDim le 2 and nDim2 le 2 then begin valid_data_or = image_or[where(~finite(image_or, /nan))] val

2020-06-27 17:26:11 806

转载 ENVI图像自动配准工具(转载)

ENVI 5.0增加了Image Registration Workflow，它是一个全新的影像配准工作流，具有自动、准确、快速的特点。它将之前版本中繁杂的参数设置步骤集成到统一的面板中，并且增加了生成种子点的影像匹配参数设置项、Harris角点算子、匹配粗差踢除算子。在少量或者无需人工干预的情况下该影像配准工作流能快速而准确的实现影像间的自动配准。功能详细操作过程启动ENVI5，在ToolBox中选择GeometricCorrection->Regi

2020-06-02 18:09:47 1777 1

原创 EndNote导入本地CAJ文件

新建一种参考文献。手动补全文献信息。右键 file attachments - attach file 选择要导入的本地CAJ文献。

2020-05-27 18:09:34 16164 1

原创 cmd调用Python脚本

cmd调用首先win + R 打开cmd窗口,cd到脚本所在的文件夹cd /d E:\pytest\银行系统练习然后输入python banksystem.py banksystem为脚本名IDL中打开cmd调用使用 & 可以一行执行多个命令spawn, 'cd /d E:\pytest\银行系统练习 & python banksystem.py'...

2020-05-21 21:24:36 450

原创 Python if name == main的简单理解

通俗来讲，if __name__ == '__main__'就是实现在程序/模块之间互相调用时声明一个运行界限的功能。程序/模块之间互相调用时if __name__ == '__main__'之前的语句将被执行，而之后的将被忽略。举个小例子：新建一个py文件，命名为’ t ‘.a = 1b = 2if __name__ == '__main__': print('a + b = {c}'.format(c=a+b))直接执行该程序得到如下输出：a + b = 3另外新建一个py文

2020-05-20 00:27:35 745

原创 MODIS产品质量控制文件使用方法

官方关于产品质量控制的说明(机翻)质量指标在生产过程中生成的CoreMetadata.0全局属性QA 中的元数据对象以及质量控制（QC）SDS中给出，或者在数据产品的产品后科学和质量检查中给出。CoreMetadata.0中的 QA元数据对象全局属性是AutomaticQualityFlag和ScienceQualityFlag及其相应的说明。根据运行LST算法期间遇到的数据条件，根据规则设置AutomaticQualityFlag。此质量检查标志的设置是完全自动化的。设置它的规则是自由的；几乎所有数

2020-05-18 00:51:13 6993 14

IDL_栅格批量转tif格式带有投影信息

空空如也