python模拟gps定位 python修改gps定位

  首先,需要获得动物跟踪研究的数据：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))

  公共区域示意图：