dtcms模板 VR全景图片 editor Java Spring forms graphics bootstrap时间轴 cpm计算 git登录命令 css鼠标悬浮样式 kubernetes安装 二分查找python pythonset python支持中文 python语言编程 python免费教程 python怎么下载 java多态 java程序实例 java怎么输出数组 php项目实例 迷宫解锁 rewritebase 快捷精灵 js倒计时 图片轮播代码 wegame更新失败 黑客攻防技术宝典 幽灵推 mysql数据库恢复 苹果内存怎么看 熊猫关键词 curdate dnf不知火刷图加点 苹果设备管理 无限弹窗代码 cgi备份还原 家长控制软件 js返回上一个页面 未来换肤
当前位置: 首页 > 学习教程  > python

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

2021/2/7 11:09:28 文章标签: 测试文章如有侵权请发送至邮箱809451989@qq.com投诉后文章立即删除

首先,需要获得动物跟踪研究的数据:Movebank   在网站中获取加拉帕戈斯信天翁的GPS定位数据,数据格式为 .csv,需要将其转换为shapefile文件,再操作数据。 数据信息: 通过location-long和location-lat字段获得x和y坐…

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

数据信息:

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

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

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

from osgeo import ogr, osr

csv_fn = r"E:\Google chrome\Download\gis with python\osgeopy data\Galapagos\Galapagos Albatrosses.csv"
shp_fn = r"E:\Google chrome\Download\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:\Google chrome\Download\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:\Google chrome\Download\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:\Google chrome\Download\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:\Google chrome\Download\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:\Google chrome\Download\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:\Google chrome\Download\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))

  公共区域示意图:

在这里插入图片描述


本文链接: http://www.dtmao.cc/news_show_2000254.shtml

附件下载

相关教程

    暂无相关的数据...

共有条评论 网友评论

验证码: 看不清楚?