Python地理数据处理 十：动物GPS点位跟踪

最新推荐文章于 2022-11-25 22:43:57 发布
最新推荐文章于 2022-11-25 22:43:57 发布
阅读量1.1k 收藏 9
点赞数 3
分类专栏： GIS with Python 文章标签： python gis
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/amyniez/article/details/113701933
版权

  首先,需要获得动物跟踪研究的数据：Movebank
在这里插入图片描述
  在网站中获取加拉帕戈斯信天翁的GPS定位数据，数据格式为 .csv，需要将其转换为shapefile文件，再操作数据。

数据信息：

在这里插入图片描述
在这里插入图片描述

  通过location-long和location-lat字段获得x和y坐标来创建一个点，并将单个本地标识符和时间戳列为属性复制。shapefile文件不能真正支持日期时间字段，所以需要将时间戳信息用字符串储存。

  1. 从一个csv文件中新建shapefile文件：

from osgeo import ogr, osr

csv_fn = r"E:\\gis with python\osgeopy data\Galapagos\Galapagos Albatrosses.csv"
shp_fn = r"E:\gis with python\osgeopy data\Galapagos\albatross_dd.shp"
sr = osr.SpatialReference(osr.SRS_WKT_WGS84_LAT_LONG)

# 创建具有两个属性字段的shapefile文件
shp_ds = ogr.GetDriverByName('ESRI Shapefile').CreateDataSource(shp_fn)
shp_lyr = shp_ds.CreateLayer('albatross_dd', sr, ogr.wkbPoint)
shp_lyr.CreateField(ogr.FieldDefn('tag_id', ogr.OFTString))
shp_lyr.CreateField(ogr.FieldDefn('timestamp', ogr.OFTString))
shp_row = ogr.Feature(shp_lyr.GetLayerDefn())

# 打开csv，循环每一行
csv_ds = ogr.Open(csv_fn)
csv_lyr = csv_ds.GetLayer()
for csv_row in csv_lyr:

    # 从csv中获取x、y坐标，并创建一个点几何图形
    x = csv_row.GetFieldAsDouble('location-long')
    y = csv_row.GetFieldAsDouble('location-lat')
    shp_pt = ogr.Geometry(ogr.wkbPoint)
    shp_pt.AddPoint(x, y)

    # 从csv中获取属性数据
    tag_id = csv_row.GetField('individual-local-identifier')
    timestamp = csv_row.GetField('timestamp')

    # 向shapefile添加数据
    shp_row.SetGeometry(shp_pt)
    shp_row.SetField('tag_id', tag_id)
    shp_row.SetField('timestamp', timestamp)
    shp_lyr.CreateFeature(shp_row)

del csv_ds, shp_ds

在这里插入图片描述
  用ArcGIS打开，显示（在非洲附近存在一个坏点，是用于数据采集过程中造成的，因此经纬度被设置为0）：

在这里插入图片描述
  2. 去除坏点（0，0）：

# 设置空间过滤条件清除坏点（0，0）
shp_ds = ogr.Open(shp_fn, True)
shp_lyr = shp_ds.GetLayer()
shp_lyr.SetSpatialFilterRect(-1, -1, 1, 1)
for shp_row in shp_lyr:
      shp_lyr.DeleteFeature(shp_row.GetFID())
shp_lyr.SetSpatialFilter(None)
shp_ds.ExecuteSQL('REPACK' + shp_lyr.GetName()) # REPACK方法永久删除点
shp_ds.ExecuteSQL('RECOMPUTE EXTENT ON' + shp_lyr.GetName()) # RECOMPUTE EXTENT ON重新计算shapefile 文件空间范围
del shp_ds

去除结果显示：

在这里插入图片描述
  3. 运用 ogr2ogr 命令行实现程序进行坐标系转换，注意：要在终端窗口中运行，并与albatross_dd.shp文件放在同一个文件夹中。
  使用兰伯特等角圆锥投影，使用它专门用于南美洲的变量。其中，o-t_srs （在输出上重新投影/转换为此SRS）和 no_def 之间的部分用来定义坐标系：

ogr2ogr -f "ESRI Shapefile" -t_srs "+proj=lcc +lat_1=-5 +lat_2=-42 +lat_0=-32 +lon_0=-60 +x_0=0 +y_0=0 +ellps=aust_SA +units=m +no_defs" E:\albatross_lambert.shp E:\albatross_dd.shp

获得以米作为单位的.shp文件:

在这里插入图片描述
  4. 选择表示一只鸟的点，从属性列中获得唯一值：

# ch7funcs模块表示从属性字段中获取唯一值
def get_unique(datasource, layer_name, field_name):
    sql = 'SELECT DISTINCT {0} FROM {1}'.format(field_name, layer_name)
    lyr = datasource.ExecuteSQL(sql)
    values = []
    for row in lyr:
        values.append(row.GetField(field_name))
    datasource.ReleaseResultSet(lyr)
    return values

  5. 计算相邻点之间的距离：

# 计算相邻点之间的距离
from osgeo import ogr
import ch7funcs

# 打开图层
# 添加一个距离字段
ds = ogr.Open(r'E:\gis with python\osgeopy data\Galapagos', True)
lyr = ds.GetLayerByName('albatross_lambert')
lyr.CreateField(ogr.FieldDefn('distance', ogr.OFTReal))

# 获取唯一的标签
# 使用ch7funcs模块
tag_ids = ch7funcs.get_unique(ds, 'albatross_lambert', 'tag_id')

# 循环遍历IDs
for tag_id in tag_ids:
    print('Processing ' + tag_id)

    # 将GPS指向限制在当前标签ID的位置
    lyr.SetAttributeFilter(
        "tag_id ='{}'".format(tag_id))
    # 保存第一次结果和点
    row = lyr.GetNextFeature()
    #row = next(lyr)
    previous_pt = row.geometry().Clone()
    previous_time = row.GetField('timestamp')

    # 循环当前标记的其余位置
    for row in lyr:
        current_time = row.GetField('timestamp')
        if current_time < previous_time:
            raise Exception('Timestamps out of order')

        # 计算到前一个点的距离并保存
        current_pt = row.geometry().Clone()
        distance = current_pt.Distance(previous_pt)
        row.SetField('distance', distance)
        lyr.SetFeature(row)

		# 保存本次计算的结果和点
        # 记住当前点，以便在处理下一个位置时，可以用作“前一个”点
        previous_pt = current_pt
        previous_time = current_time
        
del ds

  6. 找出在GPS定位点之间飞行距离最大的鸟类：

import os
import ch7funcs
from osgeo import ogr

data_dir = r"E:\gis with python\osgeopy data"
ds = ogr.Open(os.path.join(data_dir, 'Galapagos'))
tag_ids = ch7funcs.get_unique(ds, 'albatross_lambert', 'tag_id')
for tag_id in get_unique(ds, 'albatross_lambert', 'tag_id'):
    sql = """SELECT MAX(distance) FROM albatross_lambert
             WHERE tag_id = '{0}'""".format(tag_id)
    lyr = ds.ExecuteSQL(sql)
    for row in lyr:
        print('{0}: {1}'.format(tag_id, row.GetField(0)))

  运行结果：

在这里插入图片描述
  7. 查找地点间最大的飞行速度和消耗时间：

from datetime import datetime

date_format = '%Y-%m-%d %H:%M:%S.%f'
ds = ogr.Open(r'E:\gis with python\osgeopy data\Galapagos')
lyr = ds.GetLayerByName('albatross_lambert')

# 循环遍历每个标记
# 将max_speed初始化为0，并将GPS位置限制为该标记
for tag_id in ch7funcs.get_unique(ds, 'albatross_lambert', 'tag_id'):
    max_speed = 0
    lyr.SetAttributeFilter("tag_id ='{}'".format(tag_id))

    # 获取第一个点的时间戳并将其转换为datetime
    row = lyr.GetNextFeature()
    #row = next(lyr)
    ts = row.GetField('timestamp')
    previous_time = datetime.strptime(ts, date_format)

    # 循环遍历当前标记的其余位置
    for row in lyr:

        # 获取当前时间戳，转换为datetime
        # 并计算从上一个位置开始的时间
        ts = row.GetField('timestamp')
        current_time = datetime.strptime(ts, date_format)
        elapsed_time = current_time - previous_time
        hours = elapsed_time.total_seconds() / 3600

        # 计算速度
        distance = row.GetField('distance')
        speed = distance / hours

        # 保持最大速度
        max_speed = max(max_speed, speed)

    # 当循环遍历这个标签的位置时，打印出最大速度
    print('Max speed for {0}: {1}'.format(tag_id, max_speed))

  8. 用凸包多边形表示每只鸟的活动范围：

from osgeo import ogr
import ch7funcs

ds = ogr.Open(r'E:\gis with python\osgeopy data\Galapagos', True)
pt_lyr = ds.GetLayerByName('albatross_lambert')
poly_lyr = ds.CreateLayer(
    'albatross_ranges', pt_lyr.GetSpatialRef(), ogr.wkbPolygon)
id_field = ogr.FieldDefn('tag_id', ogr.OFTString)

# 创建一个足够大的面积字段来存储数据
area_field = ogr.FieldDefn('area', ogr.OFTReal)    # A
area_field.SetWidth(30)                                           
area_field.SetPrecision(4)                                  
poly_lyr.CreateFields([id_field, area_field])
poly_row = ogr.Feature(poly_lyr.GetLayerDefn())

for tag_id in ch7funcs.get_unique(ds, 'albatross_lambert', 'tag_id'):
    print('Processing ' + tag_id)
    pt_lyr.SetAttributeFilter("tag_id = '{}'".format(tag_id))

	# 用来表示位置的多点几何对象
    locations = ogr.Geometry(ogr.wkbMultiPoint)    # B
    for pt_row in pt_lyr:                                                
        locations.AddGeometry(pt_row.geometry().Clone())                 

	# 创建凸包多边形
    homerange = locations.ConvexHull()    # C
    poly_row.SetGeometry(homerange)
    poly_row.SetField('tag_id', tag_id)
    poly_row.SetField('area', homerange.GetArea())
    poly_lyr.CreateFeature(poly_row)

del ds

  创建文件：

在这里插入图片描述
  效果显示：
在这里插入图片描述
在这里插入图片描述
表示每只鸟的活动范围，显示其边界轮廓。

  9. 获取位于陆地100千米范围内的点位，创建地理区域分割的凸包多边形：

from osgeo import ogr
import ch7funcs

# 创建一个新的shp文件
ds = ogr.Open(r'E:\gis with python\osgeopy data\Galapagos', True)
pt_lyr = ds.GetLayerByName('albatross_lambert')
poly_lyr = ds.CreateLayer(
    'albatross_ranges2', pt_lyr.GetSpatialRef(), ogr.wkbPolygon)
id_field = ogr.FieldDefn('tag_id', ogr.OFTString)
area_field = ogr.FieldDefn('area', ogr.OFTReal)                          
area_field.SetWidth(30)                                                  
area_field.SetPrecision(4)
location_field = ogr.FieldDefn('location', ogr.OFTString)                                              
poly_lyr.CreateFields([id_field, area_field, location_field])
poly_row = ogr.Feature(poly_lyr.GetLayerDefn())

# 缓冲陆地多边形
land_lyr = ds.GetLayerByName('land_lambert')                             
land_row = land_lyr.GetNextFeature()                                             
land_poly = land_row.geometry().Buffer(100000)                           

for tag_id in ch7funcs.get_unique(ds, 'albatross_lambert', 'tag_id'):
    print('Processing ' + tag_id)
    pt_lyr.SetAttributeFilter("tag_id = '{}'".format(tag_id))
    pt_locations = ogr.Geometry(ogr.wkbMultiPoint)
    last_location = None
    for pt_row in pt_lyr:
        pt = pt_row.geometry().Clone()
        # 跳过不在陆地缓冲区的点
        if not land_poly.Contains(pt):                                   
            continue  
        # 弄清楚在海洋的哪一侧                                                   
        if pt.GetX() < -2800000:                                         
            location = 'island'                                          
        else:                                                            
            location = 'mainland' 
        # 如果面积改变，记录点位                                       			  
        if location != last_location:
            if pt_locations.GetGeometryCount() > 2:
                homerange = pt_locations.ConvexHull()
                poly_row.SetGeometry(homerange)
                poly_row.SetField('tag_id', tag_id)
                poly_row.SetField('area', homerange.GetArea())
                poly_row.SetField('location', last_location)
                poly_lyr.CreateFeature(poly_row)
            pt_locations = ogr.Geometry(ogr.wkbMultiPoint)
            last_location = location
        pt_locations.AddGeometry(pt)

	# 保存最后一组点
    if pt_locations.GetGeometryCount() > 2:
        homerange = pt_locations.ConvexHull()
        poly_row.SetGeometry(homerange)
        poly_row.SetField('tag_id', tag_id)
        poly_row.SetField('area', homerange.GetArea())
        poly_row.SetField('location', last_location)
        poly_lyr.CreateFeature(poly_row)

结果显示：

在这里插入图片描述
  10. 探究同一只鸟儿不同次访问岛屿之间的公共活动区域，需要计算所有岛屿访问时，公共区域所占总面积的比例：

from osgeo import ogr

ds = ogr.Open(r'E:\gis with python\osgeopy data\Galapagos')
lyr = ds.GetLayerByName('albatross_ranges2')
lyr.SetAttributeFilter("tag_id = '1163-1163' and location = 'island'")
row = lyr.GetNextFeature()
all_areas = row.geometry().Clone()
common_areas = row.geometry().Clone()
for row in lyr:
    all_areas = all_areas.Union(row.geometry())
    common_areas = common_areas.Intersection(row.geometry())
percent = common_areas.GetArea() / all_areas.GetArea() * 100
print('Percent of all area used in every visit: {0}'. format(percent))

  公共区域示意图：

在这里插入图片描述

Jackson的生态模型
关注
  • 3
    点赞
  • 踩
  • 9
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 3
    评论
专栏目录
【多点构成的面积估算】凸包算法python 实现 + 经纬度坐标 python 实现
vectorX的博客
04-09 4076
做用户聚类，有一个feature涉及到面积计算，温故一下之前了解到的凸包算法的流程，实际上，这篇文章基于的原文章，达到的精度，满足我的需求。具体实现上代码参考，因为坐标系不是笛卡尔坐标系 先上笔者参考已有的实现稍作修改后作为静态工具类使用的代码 # the class that seals the algorithm that uses convex hull to calculate ...
Scipy信号分析处理（基线漂移、滤波）（笔记01）
影子
05-15 7227
信号处理的知识： scipy.signal.detrend(): 从信号中删除线性趋势： from scipy import signal import matplotlib.pyplot as pl t = np.linspace(0, 5, 100) x = t + np.random.normal(size=100) pl.plot(t, x, linewidth=3) pl.plot(t, signal.detrend(x), linewidth=3) 效果展示： scipy.
对野生动物的GPS信息，进行数据库查询
shunrai的博客
01-01 1989
提示：文章写完后，目录可以自动生成，如何生成可参考右边的帮助文档 文章目录 前言 一、pandas是什么？ 二、使用步骤 1.引入库 2.读入数据 总结 前言 提示：这里可以添加本文要记录的大概内容： 例如：随着人工智能的不断发展，机器学习这门技术也越来越重要，很多人都开启了学习机器学习，本文就介绍了机器学习的基础内容。 提示：以下是本篇文章正文内容，下面案例可供参考 一、pandas是什么？ 示例：pandas 是基于NumPy 的一种工具，该工具是为了解决数据分
movebank-map-demo:用于从Movebank提取动物运动数据并将其显示在外部网站上的示例代码
05-06
Movebank地图演示 Movebank提供了示例代码，可用于从Movebank提取数据并使用Google Maps将其显示在外部网站上。 以下是有关其工作原理的一些详细信息： 要正常工作，您的网站必须允许您使用JavaScript代码。 某些网站提供商（例如wordpress.com）不允许这样做。 该代码有两个选项： 选项1. HTML：这里可以将整个代码发布在html页面上。 为此，必须设置Movebank中研究的权限，以便公众可以下载它们，并且必须取消选中“提示用户接受许可条款”选项。 这是因为当有人打开您的网站以查看地图时，Movebank必须将数据发送到您的网站，这与下载您的数据的网站基本相同。 如果在代码中包含用户名和密码以访问受保护的数据，则用户名和密码将在Internet上安全传递。 选项2：PHP：该代码与第一个选项非常相似，但是包含可以存储在本地服务器上PH
Python 制作国际空间站实时跟踪器
Oner.wv的专栏
10-10 598
Open Notify是一个开源项目，旨在为NASA的一些出色数据提供简单的编程接口。open-notify.org 的作者做了一些工作，以获取原始数据并将其转换为与太空和航天器有关的A...
深入实践Python-GDAL地理数据处理源码总结
03-25
本项目致力于深入实践Python-GDAL库在地理数据处理方面的应用，并对此进行了全面的源码总结。项目主要采用Python语言编写，包含总计126个文件，具体文件类型分布如下： - Python脚本文件（.py）：32个 - 图像文件...
bhulan:用于GPS数据处理的python工具箱
05-22
布兰用于GPS数据处理的开源python库可以快速处理原始gps数据，以识别给定gps轨迹内车辆运动的属性。 使用该API，您可以识别车辆的路线，停靠点，停靠持续时间，时间表和其他服务属性。基本要求 设置使用init.py文件...
基于Python的地理信息数据处理.rar
最新发布
05-07
在IT行业中，地理信息系统的应用越来越广泛，而Python作为一门强大的编程语言，因其简洁的语法和丰富的库支持，成为处理地理信息数据的理想选择。本压缩包"基于Python的地理信息数据处理.rar"提供了关于如何利用...
基于python的地理空间数据批量处理方法研究
06-14
Python语言与Arc GIS地理处理相结合,成为实现地理空间数据的批量处理快速有效的方法。文中说明了Python语言在地理空间数据处理的特点,并例举了批量定义投影、转换投影、数据裁切等空间数据批量处理的应用实例。
Python数据处理课程设计-房屋价格预测
06-14
在本“Python数据处理课程设计-房屋价格预测”项目中，我们将探讨如何运用Python编程语言和数据处理技术来预测房屋价格。房屋价格预测是经济学、统计学和机器学习领域的一个重要课题，它可以帮助政府制定政策、房...
bird-tracking::satellite:鸟类追踪-大型鸟类的GPS追踪网络
04-29
鸟类追踪 该存储库包含用于标准化LifeWatch 数据的功能，因此可以在上发布这些。 工作流程 UvA-BiTS 参考数据： 数据库→用于查询（在定义）的→可以上载到Movebank的生成的文件 GPS数据：与参考数据相同 加速度数据：与参考数据相同 奥尼特拉 参考数据：电子表格→用于查询（定义为dplyr::mutate() ）的数据的→可以上载到Movebank的生成的文件 GPS数据： 数据库→通过将选择的标签与研究相关联，实时输入Movebank 加速度数据：尚未定义 数据集 Movebank和Zenodo 标题 Movebank Zenodo 地位 BOP_RODENT-法兰德斯（比利时）啮齿类动物的猛禽（马戏团，阿西欧，布特奥） 1278021460 积极的 DELTATRACK-鲱鸥（Larus argentatus，Laridae）和小黑背鸥（Larus
GPS卫星定位—python实现
weixin_64771901的博客
11-25 4744
基于导航文件实现GPS卫星定位——python
Python实现将坐标点绘制到地图上
十八与她的博客
11-10 3971
有时我们需要直观地看到坐标点位在全球的分布情况，因此有将坐标点绘制到地图上的需求，我们可以用python实现这个功能。
【python+GIS】python实现多个栅格值提取到点，写入EXCEL表
kkkyyyxxx的博客
10-02 3207
想要将所有的温度数据（tif）的值提取到点上。
python记录点
e61015656pub的博客
12-28 96
文件编码 Unicode使用最少2个字节（1个字节=1BYTE=8bit=一个长度为8的二进制数） 来表示字母和符号等,有时候是4个字节。 UTF-8是对Unicode编码的压缩和优化，最大的特点是它采用了变长的编码方式，他不再是最少使用2个字节，而是将所有的字符进行分类。ascii码中的内容用1个字节保存、欧洲的字符用2个字节保存，东亚的字符用3个字节保存… 不同编码之间的关系...
Python使用turtle记录绘制图形经过的点的坐标-turtle.begin_poly()；turtle.end_poly()；turtle.get_poly()
liujingwei8610的专栏
08-19 1136
Python使用turtle记录绘制图形经过的点的坐标-turtle.begin_poly()；turtle.end_poly()；turtle.get_poly()
【第17届智能汽车竞赛】极速越野组——处理GPS点位的一种方法（Python-matplotlib实现手动鼠标移动坐标点的应用）
def__init__1923的博客
10-08 5038
17智能车极速越野组：GPS点位的修改是一个比较麻烦的过程，需要来回采集点位和修改程序，给调试工作带来了一定的工作量。下面这种方法可以在一定程度上方便GPS参数的修改与调试，先说一下实现的效果，可以绘制出GPS轨迹图、一键输出GPS点位程序、拖拽修改GPS点位、修改不同航段点的数量，主要是基于Python和第三方库设计出的一种方案，可以参考一下，用到的其实都是Python里面一些比较基础的知识，希望对大家能有所帮助。
python实时定位_python实现地理位置的聚类
weixin_39634237的博客
11-27 450
在地理位信息中点与点之间的位置关系是在处理数据过程中重点关注的内容,聚类算法可以将地理位置间关系直观在图形上表示出来:1、构造地理位置间距函数：defget_distance(array_1,array_2):lon_a = array_1[0]lat_a = array_1[1]lon_b = array_2[0]lat_b = array_2[1]radlat1 = radians(lat...
有关shapefile存储时间字段的问题
慕名ArcGIS
03-04 5084
很多ArcGIS的用户会很疑惑，为什么excel或者其他数据库中的时间数据数据，例如15:40:30等这样的数据，转成shapefile数据后，这个日期字段中的值会变成12:00:00？ 我们知道shapefile是基于dbase来存储的，通过dbf文件来存储属性字段，dbf支持的日期时间字段，实际上不是真正意义上的DateTime型，而是只能存储8位的日期型字段，也就是：YYY
Python空间数据处理实战：高效工具与库解析
在处理海量空间数据时，可以利用Python的并行计算能力，结合GDAL/OGR、Shapely和Rtree等库来分块处理数据，或者使用Pandas和Dask等数据处理框架进行大数据的分布式计算。通过这样的方式，可以有效地管理内存，避免一...
评论 3
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

Jackson的生态模型

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付：¥1
获取中
扫码支付

您的余额不足，请更换扫码支付或充值

打赏作者

实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值