pygplates专栏——Sample code——Create features（新建/查询特征）

Create features（新建特征）

这个例子展示了如何查询各种类型的特性。

新建常见的特征类型

这个例子展示了如何创建各种类型的特性。
这与导入类似，只是这里我们更感兴趣的是创建不同类型的特性，而不是创建通用的未分类特性类型(pygplates.FeatureType.gpml_unclassified_feature)。

由现今的几何图形创建一个海岸线特征

在这个例子中，我们从现在的几何图形中创建一个海岸线，并将其保存到一个文件中。

示例代码

import pygplates

# 首先新建一个空的海岸线特征列表
coastline_features = []
# 创建一个折线几何图形，使用部分东巴拿马海岸线的纬度/经度坐标点
eastern_panama_present_day_coastline_geometry = pygplates.PolylineOnSphere(
    [
        (-0.635538, -80.402139),
        (-0.392500, -80.500444),
        ( 0.051750, -80.079222),
        ( 0.370111, -80.049944),
        ( 0.118972, -79.998333),
        ( 0.170306, -79.948361),
        ( 0.763083, -80.099222),
        ( 0.981833, -79.650750),
        ( 0.900000, -79.649472),
        ( 1.068417, -79.430556),
        ( 1.080250, -79.177028),
        ( 1.216111, -79.061472),
        ( 1.041000, -79.014306),
        ( 1.063750, -78.928500),
        ( 1.082389, -78.982583),
        ( 1.144583, -78.897694),
        ( 1.222414, -78.932823)
    ]
)
# 根据海岸线几何图形，名称，有效时段和板块ID新建一个coastline feature
eastern_panama_coastline_feature = pygplates.Feature.create_reconstructable_feature(
    pygplates.FeatureType.gpml_coastline,
    eastern_panama_present_day_coastline_geometry,
    name="Eastern Panama, Central America",
    valid_time=(600, pygplates.GeoTimeInstant.create_distant_future()),
    reconstruction_plate_id=201
)

coastline_features.append(eastern_panama_coastline_feature)
# 保存
coastline_feature_collection = pygplates.FeatureCollection(coastline_features)
coastline_feature_collection.write("4-output_panama_coastline.gpml")

详解

需要使用一系列（比如列表）的（纬度，经度）元组来构造pygplates.PolylineOnSphere对象。当折线构造函数接收到一个2元组序列时，它将它们解释为构成折线的点的(纬度，经度)坐标。这条特殊的折线代表了今天(0Ma)巴拿马东部海岸线的部分位置。

eastern_panama_present_day_coastline_geometry = pygplates.PolylineOnSphere(
    [
        (-0.635538, -80.402139),
        (-0.392500, -80.500444),
        ( 0.051750, -80.079222),
        ( 0.370111, -80.049944),
        ( 0.118972, -79.998333),
        ( 0.170306, -79.948361),
        ( 0.763083, -80.099222),
        ( 0.981833, -79.650750),
        ( 0.900000, -79.649472),
        ( 1.068417, -79.430556),
        ( 1.080250, -79.177028),
        ( 1.216111, -79.061472),
        ( 1.041000, -79.014306),
        ( 1.063750, -78.928500),
        ( 1.082389, -78.982583),
        ( 1.144583, -78.897694),
        ( 1.222414, -78.932823)
    ]
)

使用pygplates.Feature.create_reconstructable_feature()函数创建了一个coastline feature（pygplates.FeatureType.gpml.coastline类型）
我们还传递给pygpaltes.Feature.create_reconstructable_feature函数一个现今的几何图形，一个名称，一个有效时段和一个重建板块ID。有效时段截至distant future（遥远的未来）。

eastern_panama_coastline_feature = pygplates.Feature.create_reconstructable_feature(
    pygplates.FeatureType.gpml_coastline,
    eastern_panama_present_day_coastline_geometry,
    name="Eastern Panama, Central America",
    valid_time=(600, pygplates.GeoTimeInstant.create_distant_future()),
    reconstruction_plate_id=201
)

最后保存

coastline_feature_collection = pygplates.FeatureCollection(coastline_features)
coastline_feature_collection.write('coastlines.gpml')`

备用方法

先创建一个pygplates.FeatureType.gpml_coastline空对象，然后逐一设置属性

eastern_panama_coastline_feature = pygplates.Feature(pygplates.FeatureType.gpml_coastline)
eastern_panama_coastline_feature.set_geometry(eastern_panama_present_day_coastline_geometry)
eastern_panama_coastline_feature.set_name("Eastern Panama, Central America")
eastern_panama_coastline_feature.set_valid_time(600, pygplates.GeoTimeInstant.create_distant_future())
eastern_panama_coastline_feature.set_reconstruction_plate_id(201)

进阶

pygplates.Feature.create_reconstructable_feature()函数将创建一个feature，参数type指定reconstructable features中任意类型。
如果使用了reconstructable features中任意类型的feature，那么该对象将支持gml:name，gml:description，gml:validTime和gml:reconstructionPlatedId。
在reconstructable features中有很多feature types。通过gpml:ReconstructableFeature的Inherited by features部分查看。比如说其一是gpml:TangibleFeature，其Inherited by features之一包括gpml:Coastline。这也就意味着gpml:Coastline继承自gpml:ReconsructableFeature。也就意味着它们支持相同的基础属性

从过去地质时期的几何图形中创建一个等时线特征

在这个例子中，我们从几何学中创建一个等时线，表示它在过去地质时间(而不是现在)的位置。

示例代码

import pygplates

# 加载板块运动模型
rotation_model = pygplates.RotationModel("Muller2019-Young2019-Cao2020_CombinedRotations.rot")
# 为40.1Ma的等时线坐标点创建一个折线几何图形
isochron_time_of_appearance = 40.1
isochron_geometry_at_time_of_appearance = pygplates.PolylineOnSphere(
    [
        (-57.635356,  0.765764),
        (-57.162269, -1.953176),
        (-57.916700, -2.522021),
        (-57.658576, -3.936703),
        (-58.639846, -4.849338),
        (-58.404889, -6.060713),
        (-59.390700, -6.877544),
        (-59.048499, -8.573530)
    ]
)
# 新建一个等时线feature
isochron_feature = pygplates.Feature.create_reconstructable_feature(
    pygplates.FeatureType.gpml_isochron,
    isochron_geometry_at_time_of_appearance,
    name="SOUTH AMERICAN ANTARCTIC RIDGE, SOUTH AMERICA-ANTARCTICA ANOMALY 18 IS",
    valid_time=(isochron_time_of_appearance, pygplates.GeoTimeInstant.create_distant_future()),
    reconstruction_plate_id = 201,
    conjugate_plate_id = 802,
    # 指定的几何形状不是现在的，所以它需要反向重构到现在
    reverse_reconstruct = (rotation_model, isochron_time_of_appearance)
)
# 保存
isochron_feature_collection = pygplates.FeatureCollection([isochron_feature])
isochron_feature_collection.write("4-output_isochron_40.1.gpml")

详解

pygplates.PolylineOnSphere几何图形是从(纬度，经度)元组的序列(在我们的例子中是一个列表)创建的。这是可能的，因为当折线构造函数接收到一个2元组序列时，它将它们解释为构成折线的点的(纬度，经度)坐标。

isochron_geometry_at_time_of_appearance = pygplates.PolylineOnSphere(
    [
        (-57.635356,  0.765764),
        (-57.162269, -1.953176),
        (-57.916700, -2.522021),
        (-57.658576, -3.936703),
        (-58.639846, -4.849338),
        (-58.404889, -6.060713),
        (-59.390700, -6.877544),
        (-59.048499, -8.573530)
    ])

等时线几何图形并非是现今的几何图形，所以需要反向重构创建的等时线特征到今时今日（使用reverse_reconstruct参数或pygplates.reverse_reconstruct()函数），然后才能将特征重构到任意重构时间。这是因为一个特征直到它的几何形状是现在的几何形状才算完整。
使用pygplates.Feature.create_reconstructable_feature()函数新建一个等时线特征（pygplates.FeatureType.gpml_isochron类型）。
reverse_reconstruct参数需要一个元组，包括一个板块运动模型和等时线的当前时间。
我们传递给pygplates.Feature.create_reconstructable_feature函数一个彼时彼刻的几何图形，当时的年龄（板块运动模型），名称，有效时段和重建板块ID。

isochron_feature = pygplates.Feature.create_reconstructable_feature(
    pygplates.FeatureType.gpml_isochron,
    isochron_geometry_at_time_of_appearance,
    name="SOUTH AMERICAN ANTARCTIC RIDGE, SOUTH AMERICA-ANTARCTICA ANOMALY 18 IS",
    valid_time=(isochron_time_of_appearance, pygplates.GeoTimeInstant.create_distant_future()),
    reconstruction_plate_id = 201,
    conjugate_plate_id = 802,
    # 指定的几何形状不是现在的，所以它需要反向重构到现在
    reverse_reconstruct = (rotation_model, isochron_time_of_appearance)
)

将reverse_reconstruct参数替换为pygplates.reverse_reconstruct()函数

isochron_feature = pygplates.Feature.create_reconstructable_feature(
	pygplates.FeatureType.gpml_isochron,
	isochron_geometry_at_time_of_appearance,
	name="SOUTH AMERICAN ANTARCTIC RIDGE, SOUTH AMERICA-ANTARCTICA ANOMALY 18 IS"，
	valid_time=(isochron_time_of_appearance, pygplates.GeoTimeInstant.create_distant_future()),
	reconstruction_plate_id = 201,
	conjugate_plate_id=802)
pygplates.reverse_reconstruct(isochron_feature, rotation_model, isochron_time_of_appearance)

警告

pygplates.reverse_reconstruct()要在特征的属性确定后再进行。这是必要的，因为反向重构需要查看这些属性来确定如何进行反向重构。

备用方法是在设置几何图形时反向重构：

isochron_feature = pygplates.Feature(pygplates.FeatureType.gpml.isochron)
isochron_feature.set_name("SOUTH AMERICAN ANTARCTIC RIDGE, SOUTH AMERICA-ANTARCTIC ANOMALY 18 IS")
isochron_feature.set_valid_time(isochron_time_of_appearance, pygplates.GeoTimeInstant.create_distant_future())
isochron_feature.set_reconstruction_plate_id(201)
isochron_feature.set_conjugate_plate_id(802)

isochron_feature.set_geometry(
	isochron_geometry_at_time_of_appearance)
	reverse_reconstruct=(rotation_model, isochron_time_of_appearance)))

从过去地质时期的几何图形中创建一个大洋中脊特征

示例代码

import pygplates

# 加载板块运动模型
rotation_model = pygplates.RotationModel("Muller2019-Young2019-Cao2020_CombinedRotations.rot")
# 使用过去地质时期的几何图形创建大洋中脊特征
time_of_appearance = 55.9
time_of_disappearance = 48
geometry_at_time_of_appearance = pygplates.PolylineOnSphere([...])
mid_ocean_ridge_feature = pygplates.Feature.create_tectonic_section(
    pygplates.FeatureType.gpml_mid_ocean_ridge,
    geometry_at_time_of_appearance,
    name="SOUTH ATLANTIC, SOUTH AMERICA-AFRICA",
    valid_time=(time_of_appearance, time_of_disappearance),
    left_plate=201,
    right_plate=701,
    reconstruction_method="HalfStageRotationVersion2",
    reverse_reconstruct=(rotation_model, time_of_appearance)
)

从今时今日的几何图形中创建一个俯冲带

示例代码

import pygplates

present_day_geometry = pygplates.PolylineOnSphere([...])
subduction_zone_feature = pygplates.Feature.create_reconstructable_feature(
    pygplates.FeatureType.gpml_subduction_zone,
    present_day_geometry,
    name="South America trench",
    valid_time = (200, pygplates.GeoTimeInstant.create_distant_future()),
    reconstruction_plate_id=201
)
subduction_zone_feature.set_enumeration(pygplates.PropertyName.gpml_subduction_polarity, 'Right')

我们使用pygplates.Feature.set_enumeration()将俯冲极性设置为’Right’。

创建一个虚拟地磁极特征

示例代码

import pygplates

# 极点位置和平均样本位置。
pole_position = pygplates.PointOnSphere(86.3, 168.02)
average_sample_position = pygplates.PointOnSphere(-2.91, -9.59)
# 新建虚拟地磁极特征
virtual_geomagnetic_pole_feature = pygplates.Feature.create_reconstructable_feature(
    pygplates.FeatureType.gpml_virtual_geomagnetic_pole,
    pole_position,
    name="RM: -10 -  10Ma N=10 (Dp col) Lat Range: 29.2 to -78.17 (Dm col.)",
    reconstruction_plate_id = 701
)
# 设置平均样本位置，我们需要指定它的属性名，否则它默认为极点位置并覆盖它。
virtual_geomagnetic_pole_feature.set_geometry(
    average_sample_position,
    pygplates.PropertyName.gpml_average_sample_site_position
)
# 设置平均倾角/赤纬。
virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_average_inclination, 180.16
)
virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_average_declination, 13.04
)
# 设置杆位不确定度和平均年龄。
virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_pole_a95, 3.05
)
virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_average_age, 0
)
# 保存
virtual_geomagnetic_pole_feature_collection = pygplates.FeatureCollection(virtual_geomagnetic_pole_feature)
virtual_geomagnetic_pole_feature_collection.write("4-output_virtual_geomagnetic_pole.gpml")

详解

一个虚拟地磁极特征包含两个几何图形，虚拟地磁极点位置和平均样本位置

pole_position = pygplates.PointOnSphere(86.3, 168.02)
average_sample_site_position = pygplates.PointOnSphere(-2.91, -9.59)

先通过pygplates.Feature.create_reconstructable_feature()函数创建一个虚拟地磁极点特征
这里指定的几何要素是极点位置（而不是平均样本位置）。因为虚拟地磁极点特征的默认几何要素是gpml:polePosition

virtual_geomagnetic_pole_feature = pygplates.Feature.create_reconstructable_feature(
    pygplates.FeatureType.gpml_virtual_geomagnetic_pole,
    pole_position,
    name='RM:-10 -  10Ma N= 10 (Dp col.) Lat Range: 29.2 to -78.17 (Dm col.)',
    reconstruction_plate_id=701)

我们需要单独设置平均样本位置，因为它不是默认几何要素
我们还需要指定正确的名称，否则pygplates.Feature.set_geometry()会默认设置极点位置并覆写我们设置的几何要

virtual_geomagnetic_pole_feature.set_geometry(
    average_sample_site_position,
    pygplates.PropertyName.gpml_average_sample_site_position)

接着要通过pygplates.Feature.set_double()设置一些浮点数

virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_average_inclination,
    180.16)
virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_average_declination,
    13.04)
virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_pole_a95,
    3.05)
virtual_geomagnetic_pole_feature.set_double(
    pygplates.PropertyName.gpml_average_age,
    0)

创建一个运动路径

在这个例子中，我们创建了一个运动路径特征来跟踪板块随时间的运动。

示例代码

import pygplates

# 指定今时今日非洲海岸线上的两个点（纬度/经度）
seed_points = pygplates.MultiPointOnSphere(
    [
        (-19, 12.5),
        (-28, 15.7)
    ]
)
# 重建运动轨迹的时间序列-从0-90Ma，间隔1Ma
times = range(0, 91, 1)
# 创建运动轨迹特征
motion_path_feature = pygplates.Feature.create_motion_path(
    seed_points,
    times,
    valid_time = (max(times), min(times)),
    relative_plate=201,
    reconstruction_plate_id=701
)

motion_path_feature_collection = pygplates.FeatureCollection([motion_path_feature])
motion_path_feature_collection.write("4-output_motion_path.gpml")

详解

我们指定两个点的位置在非洲的海岸线(701)。
这些是点随时间的轨迹。

seed_points = pygplates.MultiPointOnSphere(
    [
        (-19, 12.5),
        (-28, 15.7)
    ])

时间采样序列决定了运动路径的精确程度——采样的密集程度。
在这里，我们以1My的间隔从0到90Ma取样。

times = range(0, 91, 1)

我们可以通过点坐标，时间序列，有效时段和相关/重建板块ID来创建运动轨迹
这里设置有效时段为整个时间序列
点的轨迹运动基于gpml:reconstructionPlateId板块和相关的gpml:relatiePlate板块

motion_path_feature = pygplates.Feature.create_motion_path(
        seed_points,
        times,
        valid_time=(max(times), min(times)),
        relative_plate=201,
        reconstruction_plate_id=701)

创建一个流水线特征

示例代码

import pygplates

# 指定两个位于现今板块201和701之间的大洋中脊上的点
seed_points = pygplates.MultiPointOnSphere(
    [
        (-35.547600, -17.873000),
        (-46.208000, -13.623000)
    ]
)
# 取样时间序列
times = range(0, 91, 1)
# 创建流水线
flowline_feature = pygplates.Feature.create_flowline(
    seed_points,
    times,
    valid_time=(max(times),min(times)),
    left_plate=201,
    right_plate=701
)
# 保存
flowline_feature_collection = pygplates.FeatureCollection([flowline_feature])
flowline_feature_collection.write("4-output_flowline.gpml")

创建一个完整重构序列(板块运动模型)特征

在本例中，我们创建了一个表示移动板块相对于固定板块的总旋转极的时间序列的总重构序列特征。

示例代码

import pygplates
import math

# 移动板块550相对于固定板块801的有限旋转极数据。
# #数据顺序为(pole_time, pole_lat, pole_lan, pole_angle, pole_description)。
pole_data_550_rel_801 = [
    (99.0, 0.72, -179.98, 50.78, 'INA-AUS Muller et.al 2000'),
    (120.4, 10.32, -177.4, 61.12, 'INA-AUS M0 Muller et.al 2000'),
    (124.0, 11.36, -177.07, 62.54, 'INA-AUS M2 Muller et.al 2000'),
    (124.7, 11.69, -176.97, 62.99, 'INA-AUS M3 Muller et.al 2000'),
    (126.7, 12.34, -176.76, 63.95, 'INA-AUS M4 Muller et.al 2000'),
    (127.7, 12.65, -176.66, 64.42, 'INA-AUS M5 Muller et.al 2000'),
    (128.2, 12.74, -176.65, 64.63, 'INA-AUS M6 Muller et.al 2000'),
    (128.4, 12.85, -176.63, 64.89, 'INA-AUS M7 Muller et.al 2000'),
    (129.0, 13.0, -176.61, 65.23, 'INA-AUS M8 Muller et.al 2000'),
    (129.5, 13.2, -176.59, 65.67, 'INA-AUS M9 Muller et.al 2000'),
    (130.2, 13.39, -176.56, 66.1, 'INA-AUS M10 Muller et.al 2000'),
    (130.9, 13.63, -176.53, 66.66, 'INA-AUS M10N Muller et.al 2000'),
    (132.1, 13.93, -176.48, 67.4, 'INA-AUS M11 Muller et.al 2000'),
    (133.4, 14.31, -176.43, 68.33, 'INA-AUS M11A Muller et.al 2000'),
    (134.0, 14.61, -176.39, 69.09, 'INA-AUS M12 Muller et.al 2000'),
    (135.0, 14.86, -176.36, 69.73, 'INA-AUS M12A Muller et.al 2000'),
    (135.3, 15.03, -176.33, 70.19, 'INA-AUS M13 Muller et.al 2000'),
    (135.9, 15.29, -176.3, 70.89, 'INA-AUS M14 Muller et.al 2000'),
    (136.2, 15.5, -176.27, 71.44, 'INA-AUS based on closure IND-ANT Muller et.al 2000'),
    (600.0, 15.5, -176.27, 71.44, 'INA-AUS')
]
# 从极点数据创建一个有限旋转时间样本列表
pole_time_samples_550_rel_801 = [
    pygplates.GpmlTimeSample(
        pygplates.GpmlFiniteRotation(
            pygplates.FiniteRotation(pygplates.PointOnSphere(lat, lon), math.radians(angle))
        ),
        time,
        description
    ) for time, lat, lon, angle, description in pole_data_550_rel_801
]
# 时间样本需要包裹成不规则的采样属性值。
total_reconstruction_pole_550_rel_801 = pygplates.GpmlIrregularSampling(pole_time_samples_550_rel_801)
# Create the total reconstruction sequence (rotation) feature.
rotation_feature_550_rel_801 = pygplates.Feature.create_total_reconstruction_sequence(
    801,
    550,
    total_reconstruction_pole_550_rel_801,
    name="INA-AUS Muller et.al 2000"
)
# 保存
rotation_feature_550_rel_801_collectioin = pygplates.FeatureCollection([rotation_feature_550_rel_801])
rotation_feature_550_rel_801_collectioin.write("4-output_rotation_feature_550_rel_801.gpml")