锐多宝
锐多宝的朋友圈
锐多宝

使用gdal均匀筛选点矢量

作用:

通过计算各点之间的欧式距离,筛选出符合目标的、均匀发布在空间中的N个数据点。

效果示意图

运行环境

python 3.10
安装:tqdm、numpy和tqdm这三个库

完整代码

import numpy as np
from osgeo import ogr, osr
from tqdm import tqdm

# 代码作用:通过计算各点之间的欧式距离,筛选出符合目标的、均匀发布在空间中的N个数据点。



# 定义需要采样的个数
n_samples = 100
input_path = r"测试数据\村点.shp"
output_path = r"测试数据\samples.shp"

# 1. 读取原始点数据
driver = ogr.GetDriverByName('ESRI Shapefile')
inds = driver.Open(input_path, 0)
layer = inds.GetLayer()

# 2. 提取点坐标和属性
coords = []
attrs = []
for feature in layer:
    geom = feature.GetGeometryRef()
    coords.append((geom.GetX(), geom.GetY()))
    attrs.append([feature.GetField(i) for i in range(feature.GetFieldCount())])
coords = np.array(coords)
attrs = np.array(attrs)

# 3. 定义距离函数
def distance(p1, p2):
    return np.sqrt(np.sum((p1 - p2)**2))

# 4. 随机选择第一个点
idx = np.random.choice(coords.shape[0], 1)
samples = coords[idx]
sample_attrs = attrs[idx]

# 5. 选择空间均衡的采样点
for _ in tqdm(range(n_samples - 1)):
    dists = np.array([np.min(np.array([distance(p, s) for s in samples])) for p in coords])
    idx = np.argmax(dists)
    samples = np.append(samples, [coords[idx]], axis=0)
    sample_attrs = np.append(sample_attrs, [attrs[idx]], axis=0)
    coords = np.delete(coords, idx, axis=0)
    attrs = np.delete(attrs, idx, axis=0)

# 6. 将采样点转为gdal几何对象
out_samples = []
for sample in samples:
    point = ogr.Geometry(ogr.wkbPoint)
    point.AddPoint(sample[0], sample[1])
    out_samples.append(point)

# 7. 创建新的矢量层并写入采样点
out_driver = ogr.GetDriverByName('ESRI Shapefile')
out_ds = out_driver.CreateDataSource(output_path)
out_layer = out_ds.CreateLayer('samples', layer.GetSpatialRef(), ogr.wkbPoint)

# 添加属性字段
for i in range(len(layer.schema)):
    field_defn = layer.schema[i]
    out_layer.CreateField(field_defn)

# 写入采样点要素
for i, sample in enumerate(out_samples):
    feature = ogr.Feature(out_layer.GetLayerDefn())
    feature.SetGeometry(sample)
    for j, attr in enumerate(sample_attrs[i]):
        feature.SetField(j, attr)
    out_layer.CreateFeature(feature)

out_layer = None
out_ds = None
10天前
锐多宝

一条预告,一个很有用的遥感手机应用内测

引言

我之前用过吉林一号网,在上个月我又报名参加了吉林一号的免费讲座。这几天得知吉林一号网要在下周发布他们的app,我也第一时间下载了这个APP预测版本。目前该APP的功能已经开发的差不多了,提前帮大家测一测,写下这篇技术博客记录一下。

这个APP能做什么?

类似奥维的户外助手

目前支持在地图上选点,并导出。

知识课堂

这个app还提供了一些免费的遥感课程,有比较全面的遥感学习课程,看完还可以申请证书,还能参加他们的免费线下培训(我参加过一次,很赞)。

行业论坛

有挺多干货知识,从开发到应用,从webgis到遥感应用,都有涉及。

亚米级卫星底图

目前免费的版本有2022年全国一张图。即使免费版本的影像质量也挺高,提供全国免费全国0.75m分辨率影像,还可以看到影像的具体拍摄时间,如果是做遥感解译,帮助很大。

遥感影像AI分析

吉林一号网APP提供了一个户外助手的功能,可以在线解译影像、目标识别、变化检测等。我体验了一下觉得很不错,大家有兴趣可以等发布后下载玩一下这些功能。比如我试了一下,选择了我所在的区域,随意勾选了几种地物类别,过了几分钟就能获取地物覆盖图:

思考

吉林一号网APP应该是目前市面上唯一一个可以在线做遥感分析手机APP(仅限我所知道)。吉林一号网APP背靠着丰富的商业遥感卫星数据,正在创新性地制作一些非常有用的功能,比如说:

(1)全国多个时间尺度的亚米级卫星底图

(2)正在试图做一个匹配户外应用场景的功能(和奥维相比还有差距),但图源自主可控的;

(3)在APP端实现遥感影像的处理(包括地物分类、目标识别和变化检测);

(4)免费的遥感入门知识课堂分享;

(5)一个遥感行业的技术交流论坛;

我在想,如果未来这个APP可以像web端查询实时卫星影像,并购买部分区域,在后台完成影像的预处理,在线调用算法进行地物分类、目标识别和变化检测等,可能会有更广阔的市场。另外,户外助手这个功能目前相较于奥维有差距,但在这个预览版本已经能实现基本的外业调查需求,未来可期。

总结一下,吉利一号网APP非常值得下载体验,因为它拥有在线亚米底图、在线遥感影像分析功能和户外调研助手等实用功能。

参考

吉林一号.https://www.jl1mall.com/

15天前
锐多宝

2023年全国村级点位矢量制作思路与分享

引言

遇到一个网友提问,是需要所在地区的村经纬度.

因此分享一下我自己在前几个月制作的2023年全国村级点位数据的制作过程和相关数据。

处理思路

要知道全国村的分布情况,那首先我们需要知道都有哪些村?可以通过统计局的网站查询。

那有了村的名字,又怎么获取村的坐标了?地理编码

因此,我们需要先获取村名字,再进行编码,就可以获取村的分布点位矢量了。

处理步骤

获取村级名称

打开统计局的网站(https://www.stats.gov.cn/sj/tjbz/tjyqhdmhcxhfdm/2023/index.html),就可以查看各个行政区域的名称、隶属关系和区划代码了:

这个网页并没有反爬手段,而且都只是简单的html标签,所以是可以批量获取的。获取后的数据还需要经过一些数据清洗工作,才能进行地理编码。

我已经把几个年份的弄好了,有需要可以看我的这个github仓库https://github.com/ruiduobao/gaode_MAP_CUN ,在CUN这个文件夹直接下载:

因为全国同名村非常多,处理后的表格应该用全名:

地理编码

有了村的名称后,我们就可以使用高德的地理编码来获取村级矢量的分布了。

首先去高德开放平台(https://lbs.amap.com/)注册开发者账号,

获取高德的key,普通用户有5000次免费额度。我找朋友搞了一个企业账户,是30万次(有条件可以搞个企业版)。

地理编码部分,我写了个代码来获取点位的经纬度,可以参考一下:

import pandas as pd
import requests
from tqdm import tqdm
from concurrent.futures import ThreadPoolExecutor, as_completed

# 高德API的配置
GAODE_API_KEY = ''
GAODE_GEOCODE_URL = 'https://restapi.amap.com/v3/geocode/geo'

OUTPUT_DIR = r'F:\ruiduobao\map.ruiduobao.com\单个csv'

CSV_PATH = r'F:\ruiduobao\map.ruiduobao.com\1.原始数据area_code_2023_cun1.csv'

# 读取CSV文件
df = pd.read_csv(CSV_PATH)

# 定义获取经纬度的函数
def get_location(row):
    full_name = row['全名']
    url = f"{GAODE_GEOCODE_URL}?address={full_name}&key={GAODE_API_KEY}"
    try:
        response = requests.get(url)
        if response.status_code == 200:
            geocode_data = response.json()
            if geocode_data and geocode_data['geocodes']:
                location = geocode_data['geocodes'][0]['location']
                return location.split(',')
            else:
                return ('获取失败', '获取失败')
        else:
            return ('获取失败', '获取失败')
    except requests.RequestException:
        return ('获取失败', '获取失败')

# 使用线程池来加速请求处理
with ThreadPoolExecutor(max_workers=30) as executor:  # 可以调整max_workers参数来增加或减少线程数
    # 将每一行数据提交给线程池
    futures = [executor.submit(get_location, row) for index, row in df.iterrows()]
    
    # 初始化进度条
    progress = tqdm(as_completed(futures), total=len(futures), desc="Processing", unit="row")
    
    # 收集结果,并且在进度条中更新进度
    for future in progress:
        result = future.result()  # 获取线程结果
        index = futures.index(future)  # 获取任务对应的索引
        try:
            longitude, latitude = result
            df.at[index, '经度'] = longitude
            df.at[index, '纬度'] = latitude
        except Exception as exc:
            df.at[index, '经度'] = '获取失败'
            df.at[index, '纬度'] = '获取失败'
            print(f"{index} generated an exception: {exc}")

# 写回更新后的CSV文件
df.to_csv(CSV_PATH, index=False)

火星坐标系转换

由于我国行业政策原因,通过地理编码获取的经纬度是火星坐标系,有偏移,下面绿色点是通过高德获取到的村级点位(有偏移),红色是wgs84坐标系。

把火星坐标转为WGS84的工具,我比较推荐Geohy,但是这个工具的官方在线版本已经不提供转wgs84

但我们可以在github上找未改动的geohey的版本(https://github.com/ruiduobao/qgis-geohey-toolbox),下载插件再本地安装,即可使用。

全国大概是60多万个村级点,等待10分钟左右即可全部转为wgs84的矢量。

数据分享

这个是我已经处理好的村级点位数据。

百度云数据链接:https://pan.baidu.com/s/1xYjoA5ZVFJAytQ110ZNQaw?pwd=4444
提取码:4444

总结

总结一下技术路线,要获取全国村级点位矢量,需要:

(1)统计局先获得全国各级行政区名称和代码名称(用爬虫或者github)

(2)用高德的地理编码,获取每一个村级点的经纬度(火星坐标系)

(3)将火星坐标系再转为WGS84坐标系,有多种工具可以选择,这里我使用的qgis插件geohey。

2024年04月10日
锐多宝

网格矢量如何计算莫兰指数

引言

遇到一个问题,计算矢量网格的莫兰指数。

概念解释

莫兰指数

莫兰指数(Moran's Index)是一种空间自相关指标,用于衡量空间数据的相似性和聚集程度。它可以用来描述一个区域与其邻近区域之间的属性值的相关性。莫兰指数的取值范围通常在-1到1之间。

  • 当莫兰指数接近1时,表示空间数据呈现出正相关,即相似的值倾向于聚集在一起。
  • 当莫兰指数接近-1时,表示空间数据呈现出负相关,即不同的值倾向于聚集在一起。
  • 当莫兰指数接近0时,表示空间数据呈现出随机分布,没有明显的空间自相关性。

knearst=4?

knearst=4矩阵是一种空间权重矩阵,用于定义空间数据中每个观测点的邻域。在这种矩阵中,每个观测点的邻域由其最近的4个点组成。

解决思路

计算矢量数据中每个要素(网格)的局部莫兰指数,并将计算结果添加到矢量数据的属性表中。我做了一个示意矢量,如图所示:

因为需要涉及到矢量数据的操作,这里我们使用gdal

还涉及到莫兰指数,我们使用pysal,这个包用于空间权重矩阵的构建、空间自相关指标的计算、空间回归模型的估计等。

初始化和读取矢量数据

import numpy as np
import pysal
from osgeo import ogr

driver = ogr.GetDriverByName('ESRI Shapefile')
SHP_PATH = r"矢量数据.shp"
dataset = driver.Open(SHP_PATH, 1) 
layer = dataset.GetLayer()
  1. 使用 ogr 库打开矢量数据文件(ESRI Shapefile),以读写模式打开。
  2. 获取矢量数据的图层。

提取属性值和坐标

values = []

coords = []
for feature in layer:

geom = feature.GetGeometryRef()
centroid = geom.Centroid()
coords.append([centroid.GetX(), centroid.GetY()])
values.append(feature.GetField('singlearea'))

values = np.array(values)
coords = np.array(coords)

  1. 遍历图层中的每个要素(feature)。
  2. 获取要素的几何体(geometry),并计算其质心坐标。
  3. 将质心坐标添加到 coords 列表中。
  4. 将指定字段('singlearea')的属性值添加到 values 列表中。
  5. 将属性值和坐标转换为 NumPy 数组。

创建权重矩阵

knn = pysal.lib.weights.KNN(coords, k=4)
knn.transform = 'r'
  1. 使用 pysal 库的 KNN 函数创建 k 最近邻权重矩阵,设置 k=4
  2. 对权重矩阵进行行标准化。

计算局部莫兰指数

local_moran = pysal.explore.esda.Moran_Local(values, knn)
print("局部莫兰指数:", local_moran.Is)

标准化局部莫兰指数

min_value = np.min(local_moran.Is)
max_value = np.max(local_moran.Is)
normalized_local_moran = (local_moran.Is - min_value) / (max_value - min_value) * 2 - 1
print("标准化后的局部莫兰指数:", normalized_local_moran)

  1. 使用 pysal 库的 Moran_Local 函数计算每个网格的局部莫兰指数。
  2. 打印计算得到的局部莫兰指数。

将局部莫兰指数添加到矢量数据属性表

lisa_field = ogr.FieldDefn('LISA_I', ogr.OFTReal)
layer.CreateField(lisa_field)

dataset = None
dataset = driver.Open(SHP_PATH, 1)
layer = dataset.GetLayer()

for i in range(layer.GetFeatureCount()):

feature = layer.GetFeature(i)
feature.SetField('LISA_I', float(local_moran.Is[i]))
layer.SetFeature(feature)</code></pre><ol><li>创建一个新的字段('LISA_I')来存储局部莫兰指数。</li><li>重新打开矢量数据集并获取图层。</li><li>遍历图层中的每个要素。</li><li>使用 <code>layer.GetFeature(i)</code> 获取要素,并将对应的局部莫兰指数赋值给新字段。</li><li>更新要素的属性表。</li></ol><h3>关闭数据集并销毁数据源</h3><pre><code>dataset.Destroy()

dataset = None
print("局部莫兰指数已成功添加到矢量数据属性表中。")

  1. 关闭矢量数据集。
  2. 销毁数据源以释放资源。
  3. 打印提示信息,表示局部莫兰指数已成功添加到矢量数据的属性表中。

完整代码

import numpy as np
import pysal
from osgeo import ogr

打开矢量数据文件(以读写模式打开)

driver = ogr.GetDriverByName('ESRI Shapefile')
SHP_PATH = r"矢量数据 - 副本.shp"
dataset = driver.Open(SHP_PATH, 1)
layer = dataset.GetLayer()

提取属性值和坐标

values = []
coords = []
for feature in layer:

geom = feature.GetGeometryRef()
centroid = geom.Centroid()
coords.append([centroid.GetX(), centroid.GetY()])
values.append(feature.GetField('cenlan'))

将属性值和坐标转换为NumPy数组

values = np.array(values)
coords = np.array(coords)

创建k最近邻权重矩阵(knearst=4)

knn = pysal.lib.weights.KNN(coords, k=4)

行标准化权重矩阵

knn.transform = 'r'

计算每个网格的局部莫兰指数

local_moran = pysal.explore.esda.Moran_Local(values, knn)
print("局部莫兰指数:", local_moran.Is)

标准化局部莫兰指数

min_value = np.min(local_moran.Is)
max_value = np.max(local_moran.Is)
normalized_local_moran = (local_moran.Is - min_value) / (max_value - min_value) * 2 - 1
print("标准化后的局部莫兰指数:", normalized_local_moran)

将标准化后的局部莫兰指数添加到矢量数据属性表,使用有效的字段名称

lisa_field = ogr.FieldDefn('LISA_I', ogr.OFTReal)
layer.CreateField(lisa_field)

重新打开数据集并获取图层

dataset = None
dataset = driver.Open(SHP_PATH, 1)
layer = dataset.GetLayer()

使用 layer.GetFeature(i) 获取要素并更新,使用更新后的字段名称

for i in range(layer.GetFeatureCount()):

feature = layer.GetFeature(i)
feature.SetField('LISA_I', float(normalized_local_moran[i]))
layer.SetFeature(feature)

关闭数据集并销毁数据源

dataset.Destroy()
dataset = None

print("标准化后的局部莫兰指数已成功添加到矢量数据属性表中。")

效果展示

运行完代码,效果为:

总结

使用gdal负责空间数据处理,使用pysal完成莫兰指数的计算,然后把计算结果写入到属性表里,

2024年04月09日
锐多宝

QGIIS制作萤火图

今天研究了一下萤火图,以为很难,实际上就两步就设置完成了,具体可以看我的截图:

首先要准备一个暗色背景,在我这里用的ESRI的暗色背景图,然后再准备一个点数据(我这里采用的是北京市的博物馆分布图)。

选中点图层,按照下面的设置进行调整:

在绘制效果中,点击进入,修改为:

最后就能生成萤火虫般的点分布图了

总结

简单的两步就能使用QGIS创建点的萤火图了,就是利用点的外发光以及模糊半径来进行显示。快来试试。

参考

https://X.com/helenmakesmaps

2024年03月31日
锐多宝 .