python 对陆地数据进行掩膜的两种方法

简介: 一般我们分析海洋或者气象相关数据时,经常会绘制散点图分析两重变量之间的相关性。近期在分析数据时,发现有个容易出问题的情况:

前提杂谈



一般我们分析海洋或者气象相关数据时,经常会绘制散点图分析两重变量之间的相关性。近期在分析数据时,发现有个容易出问题的情况:


  • 如果我们只需要分析海洋上的某一变量的变化时,最好先将陆地数据进行掩膜,防止因为地形等原因导致陆地上的数据存在误差,影响我们的分析结果。
    基于上述需求,在这里介绍两种掩膜的方法,这里选择将陆地进行掩膜,掩膜海洋的办法类似。
  • 一种是在任何windows系统上的python自带的库进行掩膜,另一种需要通过geopandas+Salem库进行掩膜。其中,后者往往库在Windows上难以安装,一般建议在linux系统进行处理。


方法



  • 1、使用python自带的库global_land_mask进行掩膜 ( Windows/Linux都适用)
  • 2、使用geopandas+Salem进行掩膜 (建议在Linux系统上使用)


方法1:使用global_land_mask进行掩膜


这里主要通过global_land_mask中的globe.is_ocean()判断数据点的所在的经纬度坐标是否在海洋。如果是陆地,就进行掩膜,反之就不处理。下面封装一个掩膜的函数,方便后续调用处理。


所需要的库主要如下所示:

from global_land_mask import globe
import cmaps
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from cartopy.mpl.ticker import LongitudeFormatter, LatitudeFormatter


封装掩膜的函数代码如下,其中,对于不同文件的经度名称也进行了考虑设置。

def mask_land(ds, label='land', lonname='lon'):
    if lonname == 'lon':
        lat = ds.lat.data
        lon = ds.lon.data
        if np.any(lon > 180):
            lon = lon - 180
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
            temp = []
            temp = mask[:, 0:(len(lon) // 2)].copy()
            mask[:, 0:(len(lon) // 2)] = mask[:, (len(lon) // 2):]
            mask[:, (len(lon) // 2):] = temp
        else:
            lons, lats = np.meshgrid(lon, lat)# Make a grid
            mask = globe.is_ocean(lats, lons)# Get whether the points are on ocean.
        ds.coords['mask'] = (('lat', 'lon'), mask)
    elif lonname == 'longitude':
        lat = ds.latitude.data
        lon = ds.longitude.data
        if np.any(lon > 180):
            lon = lon - 180
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
            temp = []
            temp = mask[:, 0:(len(lon) // 2)].copy()
            mask[:, 0:(len(lon) // 2)] = mask[:, (len(lon) // 2):]
            mask[:, (len(lon) // 2):] = temp
        else:
            lons, lats = np.meshgrid(lon, lat)
            mask = globe.is_ocean(lats, lons)
        lons, lats = np.meshgrid(lon, lat)
        mask = globe.is_ocean(lats, lons)
        ds.coords['mask'] = (('latitude', 'longitude'), mask)
    if label == 'land':
        ds = ds.where(ds.mask == True)
    elif label == 'ocean':
        ds = ds.where(ds.mask == False)
    return ds


这个函数调用起来非常简单,只需要输入你的数据,选择掩膜的区域是海洋和陆地,指定经度名称即可。


下面以月均海温资料为例,展示掩膜前后的效果:


2bd185809eb54bf9a580c8445d0a0a0a.png

可以发现,掩膜成功,效果还是不错的。基本满足我们下一步的处理需求。绘图代码如下:

infile = xr.open_dataset(r"sst.mnmean.nc")
infile = mask_land(infile,'land','lon')
infile = infile.transpose("time","lat","lon",...)
def make_map(ax, title):
    # set_extent  set crs
    ax.set_extent(box, crs=ccrs.PlateCarree())
    # land = cfeature.NaturalEarthFeature('physical',
    #                                     'land',
    #                                     scale,
    #                                     edgecolor='face')
    ax.set_xticks(np.arange(box[0], box[1] + xstep, xstep),
                  crs=ccrs.PlateCarree())  # [ )
    ax.set_yticks(np.arange(box[2], box[3] + ystep, ystep),
                  crs=ccrs.PlateCarree())
    plt.tick_params(labelsize=15)
    lon_formatter = LongitudeFormatter(zero_direction_label=False)  
    lat_formatter = LatitudeFormatter()
    ax.gridlines()
    ax.xaxis.set_major_formatter(lon_formatter)
    ax.yaxis.set_major_formatter(lat_formatter)
    ax.set_ylabel('Latitude', fontsize=15)
    ax.set_xlabel('Longitute', fontsize=15)
    ax.set_title(title, fontsize=23, loc='left')
    return ax
lat = infile['lat'].data
lon = infile['lon'].data
sst = infile['sst'].data
sst_m=np.nanmean(sst,axis=0)
sst_ori = xr.open_dataset(r"sst.mnmean.nc").sst
sst_ori = np.nanmean(sst_ori,axis=0)
box = [0, 361, -90, 85]  
scale = '50m'            
xstep, ystep = 45, 35   
# cmap=plt.get_cmap('rainbow')#'RdYlBu_r'
cmap=cmaps.NCV_jet
fig=plt.figure(figsize=(20,16))
ax=fig.add_subplot(121,projection=ccrs.PlateCarree(central_longitude=180))
make_map(ax,'SST-Mask data')
ax.contourf(lon,lat,sst_m,cmap=cmap,transform=ccrs.PlateCarree(),zorder=2)
ax2=fig.add_subplot(122,projection=ccrs.PlateCarree(central_longitude=180))
make_map(ax2,'SST-origin data')
plot=ax2.contourf(lon,lat,sst_ori,cmap=cmap,transform=ccrs.PlateCarree(),zorder=2)
ax3=fig.add_axes([0.27,0.3,0.4,0.017])
cb=fig.colorbar(plot,cax=ax3,shrink=0.9,ticks=[0,10,20,30],pad=0.04,aspect=15,orientation='horizontal')
cb.ax.tick_params(labelsize=10)
cb.ax.set_title('$°C$', fontsize=20, loc='center')


方法2:使用geopandas+salem库进行掩膜


关于salem和geopandas在Windows系统的安装可以查看以下帖子,可能有用:

salem安装

geopandas安装


使用geopandas+salem进行掩膜相比上述方法,大大减少了代码量,非常方便简洁,只需要提供数据和一个海洋的shp文件即可。缺点在于geopanda和Salem在windows系统的安装较为困难,很难成功安装使用。但是,在linux系统下则相对来说比较容易安装。这里也在Linux系统下以对垂直速度的数据的掩膜进行演示:

我这里的数据的经度排序是0~360,先将其转为-180~180排列,方便我后续处理使用:


import salem
import geopandas as geo
import xarray as xr
ds=xr.open_dataset("omiga.nc")
# change 0-360 to -180 ~ 180
lon_name = 'lon'
ds['longitude_adjusted'] = xr.where(
    ds[lon_name] > 180,
    ds[lon_name] - 360,
    ds[lon_name])
ds = (ds.swap_dims({lon_name: 'longitude_adjusted'}).sel(**{'longitude_adjusted': sorted(ds.longitude_adjusted)}).drop(lon_name))
ds = ds.rename({'longitude_adjusted': lon_name})
# mask
shp_path="shp/10m_physical/ne_10m_ocean_scale_rank.shp"
shp=geo.read_file(shp_path)
sst=ds.salem.roi(shape=shp)
sst.to_netcdf(path="omiga_mask.nc")


绘制填色图查看掩膜结果:


294c8ac69f6a4cff91e8268b8700e505.png

总体上还不错,绘图代码这里就不贴了,与方法1中的代码一致。

虽然不同的掩膜方法掩膜的结果可能不同,但是基本上满足我们后续的计算需求,只需要统一掩膜方法即可。有兴趣的xd可以尝试起来。


                        一个努力学习python的海洋人
                        水平有限,欢迎指正!!!
                          欢迎评论、收藏、点赞、转发、关注。
                          关注我不后悔,记录学习进步的过程~~


相关文章
|
1天前
|
算法 Serverless 数据处理
从集思录可转债数据探秘:Python与C++实现的移动平均算法应用
本文探讨了如何利用移动平均算法分析集思录提供的可转债数据,帮助投资者把握价格趋势。通过Python和C++两种编程语言实现简单移动平均(SMA),展示了数据处理的具体方法。Python代码借助`pandas`库轻松计算5日SMA,而C++代码则通过高效的数据处理展示了SMA的计算过程。集思录平台提供了详尽且及时的可转债数据,助力投资者结合算法与社区讨论,做出更明智的投资决策。掌握这些工具和技术,有助于在复杂多变的金融市场中挖掘更多价值。
22 12
|
6天前
|
数据采集 数据安全/隐私保护 Python
从零开始:用Python爬取网站的汽车品牌和价格数据
在现代化办公室中,工程师小李和产品经理小张讨论如何获取懂车帝网站的汽车品牌和价格数据。小李提出使用Python编写爬虫,并通过亿牛云爬虫代理避免被封禁。代码实现包括设置代理、请求头、解析网页内容、多线程爬取等步骤,确保高效且稳定地抓取数据。小张表示理解并准备按照指导操作。
从零开始:用Python爬取网站的汽车品牌和价格数据
|
11天前
|
JSON 数据可视化 API
Python 中调用 DeepSeek-R1 API的方法介绍,图文教程
本教程详细介绍了如何使用 Python 调用 DeepSeek 的 R1 大模型 API,适合编程新手。首先登录 DeepSeek 控制台获取 API Key,安装 Python 和 requests 库后,编写基础调用代码并运行。文末包含常见问题解答和更简单的可视化调用方法,建议收藏备用。 原文链接:[如何使用 Python 调用 DeepSeek-R1 API?](https://apifox.com/apiskills/how-to-call-the-deepseek-r1-api-using-python/)
|
26天前
|
数据挖掘 数据处理 开发者
Python3 自定义排序详解:方法与示例
Python的排序功能强大且灵活,主要通过`sorted()`函数和列表的`sort()`方法实现。两者均支持`key`参数自定义排序规则。本文详细介绍了基础排序、按字符串长度或元组元素排序、降序排序、多条件排序及使用`lambda`表达式和`functools.cmp_to_key`进行复杂排序。通过示例展示了如何对简单数据类型、字典、类对象及复杂数据结构(如列车信息)进行排序。掌握这些技巧可以显著提升数据处理能力,为编程提供更强大的支持。
32 10
|
28天前
|
人工智能 自然语言处理 算法
随机的暴力美学蒙特卡洛方法 | python小知识
蒙特卡洛方法是一种基于随机采样的计算算法,广泛应用于物理学、金融、工程等领域。它通过重复随机采样来解决复杂问题,尤其适用于难以用解析方法求解的情况。该方法起源于二战期间的曼哈顿计划,由斯坦尼斯拉夫·乌拉姆等人提出。核心思想是通过大量随机样本来近似真实结果,如估算π值的经典示例。蒙特卡洛树搜索(MCTS)是其高级应用,常用于游戏AI和决策优化。Python中可通过简单代码实现蒙特卡洛方法,展示其在文本生成等领域的潜力。随着计算能力提升,蒙特卡洛方法的应用范围不断扩大,成为处理不确定性和复杂系统的重要工具。
69 21
|
1月前
|
数据采集 Web App开发 数据可视化
Python用代理IP获取抖音电商达人主播数据
在当今数字化时代,电商直播成为重要的销售模式,抖音电商汇聚了众多达人主播。了解这些主播的数据对于品牌和商家至关重要。然而,直接从平台获取数据并非易事。本文介绍如何使用Python和代理IP高效抓取抖音电商达人主播的关键数据,包括主播昵称、ID、直播间链接、观看人数、点赞数和商品列表等。通过环境准备、代码实战及数据处理与可视化,最终实现定时任务自动化抓取,为企业决策提供有力支持。
|
2月前
|
数据采集 存储 XML
python实战——使用代理IP批量获取手机类电商数据
本文介绍了如何使用代理IP批量获取华为荣耀Magic7 Pro手机在电商网站的商品数据,包括名称、价格、销量和用户评价等。通过Python实现自动化采集,并存储到本地文件中。使用青果网络的代理IP服务,可以提高数据采集的安全性和效率,确保数据的多样性和准确性。文中详细描述了准备工作、API鉴权、代理授权及获取接口的过程,并提供了代码示例,帮助读者快速上手。手机数据来源为京东(item.jd.com),代理IP资源来自青果网络(qg.net)。
|
2月前
|
数据采集 Web App开发 监控
Python爬虫:爱奇艺榜单数据的实时监控
Python爬虫:爱奇艺榜单数据的实时监控
|
2月前
|
安全
Python-打印99乘法表的两种方法
本文详细介绍了两种实现99乘法表的方法:使用`while`循环和`for`循环。每种方法都包括了步骤解析、代码演示及优缺点分析。文章旨在帮助编程初学者理解和掌握循环结构的应用,内容通俗易懂,适合编程新手阅读。博主表示欢迎读者反馈,共同进步。
|
2月前
|
数据采集 分布式计算 大数据
构建高效的数据管道:使用Python进行ETL任务
在数据驱动的世界中,高效地处理和移动数据是至关重要的。本文将引导你通过一个实际的Python ETL(提取、转换、加载)项目,从概念到实现。我们将探索如何设计一个灵活且可扩展的数据管道,确保数据的准确性和完整性。无论你是数据工程师、分析师还是任何对数据处理感兴趣的人,这篇文章都将成为你工具箱中的宝贵资源。

热门文章

最新文章

推荐镜像

更多