小编stm*_*4tt的帖子

我已经将此问题作为在不规则网格问题上绘制数据的有效方法的一部分,但一般反馈是将原始问题拆分为更易于管理的块.因此,这个新问题.

我使用组织在不规则二维网格上的卫星数据,其尺寸为扫描线(沿轨道尺寸,即Y轴)和地面像素(跨越轨道尺寸,即X轴).每个中心像素的纬度和经度信息存储在辅助坐标变量中,以及四个角坐标对(纬度和经度坐标在WGS84参考椭球上给出).

让我们构建一个玩具数据集,包括12x10可能不规则的网格和相关的表面温度测量.

library(pracma) # for the meshgrid function
library(ggplot2)

num_sl <- 12 # number of scanlines
num_gp <- 10 # number of ground pixels
l <- meshgrid(seq(from=-20, to=20, length.out = num_gp), 
              seq(from=30, to=60, length.out = num_sl))

lon <- l[[1]] + l[[2]]/10
lat <- l[[2]] + l[[1]]/10

data <- matrix(seq(from = 30, to = 0, length.out = num_sl*num_gp), 
               byrow = TRUE, nrow = num_sl, ncol = num_gp) +
  matrix(runif(num_gp*num_sl)*6, nrow = num_sl, ncol = …

我试图在Python中提出一种通用的方法来识别在一组计划的航天器机动过程中发生的俯仰旋转.您可以将其视为移位检测问题的特定情况.

让我们考虑solar_elevation_angle我的测量集中的变量,确定从航天器仪器测量的太阳仰角.对于那些可能想要使用数据的人,我在这里保存了solar_elevation_angle.txt文件.

import numpy as np
import matplotlib.pyplot as plt
from matplotlib import gridspec
from scipy.signal import argrelmax
from scipy.ndimage.filters import gaussian_filter1d

solar_elevation_angle = np.loadtxt("solar_elevation_angle.txt", dtype=np.float32)

fig, ax = plt.subplots()    
ax.set_title('Solar elevation angle')
ax.set_xlabel('Scanline')
ax.set_ylabel('Solar elevation angle [deg]')
ax.plot(solar_elevation_angle)
plt.show()

扫描线是我的时间维度.斜率变化的四个点识别航天器俯仰旋转.

正如您所看到的,航天器机动区域外的太阳高度角演变与时间的关系几乎是线性的,对于这个特定的航天器而言应该始终如此(主要故障除外).

请注意,在每次航天器操作期间,虽然在我的角度值组中离散,但斜率变化显然是连续的.这意味着:对于每次操作,尝试定位已进行机动的单个扫描线并不是真的有意义.我的目标是为每次操作识别扫描线范围内的"代表性"扫描线,该扫描线限定机动发生的时间间隔(例如,中间值或左边界).

一旦我获得了一组"代表性"扫描线索引,其中所有操作都已发生,我可以使用这些索引粗略估计操纵持续时间,或自动在标绘上放置标签.

到目前为止我的解决方案是:

1. 使用计算太阳高度角的二阶导数 np.gradient.
2. 计算绝对值并剪切生成的曲线.剪切是必要的,因为我假设线性段中的离散化噪声,这将严重影响点4中"实际"局部最大值的识别.
3. 对所得曲线应用平滑,以消除多个峰.我正在使用scipy的1d高斯滤波器,其具有反复试验的sigma值.
4. 识别局部最大值.

这是我的代码:

fig = plt.figure(figsize=(8,12))
gs = gridspec.GridSpec(5, 1) 

ax0 = plt.subplot(gs[0])
ax0.set_title('Solar elevation angle')
ax0.plot(solar_elevation_angle)

solar_elevation_angle_1stdev = np.gradient(solar_elevation_angle)
ax1 = plt.subplot(gs[1]) …

我使用组织在不规则二维网格上的卫星数据,其尺寸为扫描线(沿轨道尺寸)和地面像素(跨越轨道尺寸).每个中心像素的纬度和经度信息存储在辅助坐标变量中,以及四个角坐标对(纬度和经度坐标在WGS84参考椭球上给出).数据存储在netCDF4文件中.

我想要做的是在投影地图上有效地绘制这些文件(可能还有文件的组合 - 下一步!).

我的做法,到目前为止,灵感来自杰里米的Voisey的回答这个问题,一直在打造我的兴趣可变连杆的像素边界的数据帧,并使用ggplot2与geom_polygon用于实际的情节.

让我来说明我的工作流程,并提前为天真的方法道歉:我刚开始用一周或两周的R编码.

注意

要完全重现问题:
1.下载两个数据帧:so2df.Rda(22M)和pixel_corners.Rda(26M)
2.在您的环境中加载它们,例如

so2df <- readRDS(file="so2df.Rda")
pixel_corners <- readRDS(file="pixel_corners.Rda")

3. 跳转到"合并数据帧"步骤.

初始设置

我要从我的文件中读取数据和纬度/经度边界.

library(ncdf4)
library(ggplot2)
library(ggmap) 
# set path and filename
ncpath <- "/Users/stefano/src/s5p/products/e1dataset/L2__SO2/"
ncname <- "S5P_OFFL_L2__SO2____20171128T234133_20171129T003956_00661_01_022943_00000000T000000"  
ncfname <- paste(ncpath, ncname, ".nc", sep="")
nc <- nc_open(ncfname)

# save fill value and multiplication factors
mfactor = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column", 
                    "multiplication_factor_to_convert_to_DU")
fillvalue = ncatt_get(nc, "PRODUCT/sulfurdioxide_total_vertical_column", 
                      "_FillValue")

# read the SO2 total column variable
so2tc <- …

我对Cartopy渲染Google瓷砖感到有些困惑。与标准的Google地图外观相比，该地图看起来非常差。

示例（来自https://ocefpaf.github.io/python4oceanographers/blog/2015/06/22/osm/的代码）：

import matplotlib.pyplot as plt

import cartopy.crs as ccrs
from cartopy.io import shapereader
from cartopy.mpl.gridliner import LONGITUDE_FORMATTER, LATITUDE_FORMATTER

def make_map(projection=ccrs.PlateCarree()):
    fig, ax = plt.subplots(figsize=(9, 13),
                           subplot_kw=dict(projection=projection))
    gl = ax.gridlines(draw_labels=True)
    gl.xlabels_top = gl.ylabels_right = False
    gl.xformatter = LONGITUDE_FORMATTER
    gl.yformatter = LATITUDE_FORMATTER
    return fig, ax
import cartopy.io.img_tiles as cimgt

extent = [-39, -38.25, -13.25, -12.5]

request = cimgt.GoogleTiles()

fig, ax = make_map(projection=request.crs)
ax.set_extent(extent)

ax.add_image(request, 10)

产生：

与链接的网站上显示的同一张图片相比，这看起来非常糟糕-查看文字标签和街道编号的像素化渲染：

更改缩放级别似乎并不能改善这种情况。

这是由Cartopy和googletiles（）渲染的地图上的另一个示例：

在Google地图中显示的同一张地图

有谁知道这个奇怪问题的起因是什么以及如何解决？

我使用在不规则二维网格上组织的卫星数据，其尺寸为扫描线（沿轨道尺寸）和地面像素（跨轨道尺寸）。每个地面像素的纬度和经度信息存储在辅助坐标变量中。

给定一个（纬度，经度）点，我想识别我的数据集上最接近的地面像素。

让我们构建一个 10x10 的玩具数据集：

import numpy as np
import xarray as xr
import cartopy.crs as ccrs
import matplotlib.pyplot as plt
%matplotlib inline

lon, lat = np.meshgrid(np.linspace(-20, 20, 10), 
                       np.linspace(30, 60, 10))
lon += lat/10
lat += lon/10
da = xr.DataArray(data = np.random.normal(0,1,100).reshape(10,10), 
                  dims=['scanline', 'ground_pixel'],
                  coords = {'lat': (('scanline', 'ground_pixel'), lat),
                            'lon': (('scanline', 'ground_pixel'), lon)})

ax = plt.subplot(projection=ccrs.PlateCarree())
da.plot.pcolormesh('lon', 'lat', ax=ax, cmap=plt.cm.get_cmap('Blues'), 
                   infer_intervals=True);
ax.scatter(lon, lat, transform=ccrs.PlateCarree())
ax.coastlines()
ax.gridlines(draw_labels=True)
plt.tight_layout()

请注意，纬度/经度坐标标识中心像素，像素边界由 xarray 自动推断。

现在，假设我想识别距离罗马最近的地面像素。

到目前为止，我想到的最好方法是在堆叠的扁平纬度/经度数组上使用 scipy 的 kdtree：

from scipy import …