0

我是形状文件的新手。我在太平洋工作，经常遇到与日期变更相关的问题。

目前，我能够为 从https://www.marineregions.org/gazetteer.php?p=details&id=21483下载的太平洋longhurst 形状文件绘制 Longhurst 省的一半

需要下载 longhurst shapefile 到 longhurst/

from mpl_toolkits.basemap import Basemap
import geopandas as gp
m = Basemap(llcrnrlon=120.,llcrnrlat=-40,urcrnrlon=290,urcrnrlat=40,
                resolution='l',projection='merc')
m.readshapefile('longhurst/Longhurst_world_v4_2010',name='longhurst')

这会产生这样的东西。我想为整个太平洋做这件事。不只是西方：

破碎的省份地图

然后我发现我可以用 geopandas 打开形状文件。 provinces=gp.read_file('datasets/longhurst/')

我无法弄清楚如何修改经度。我只需要+360，或者我认为更复杂的它们，并且可以绘制两个版本以获得跨越日期线的范围。这有什么智慧吗？谢谢。

在 xarray 中，以下工作可以改变经度......但不确定如何应用于这种情况。

dat= dat.assign_coords(lon=(dat.lon % 360)).roll(lon=dat.dims['lon']),roll_coords=False).sortby('lon')

4

2 回答 2

0

我已经编辑了我在另一个问题中回答的代码，以解决您的一些问题（并非所有人都可能）。它是可运行的（使用某些 shapefile）并生成附加图。我希望该代码对您进一步调查有用。

import cartopy.crs as ccrs
import geopandas
import matplotlib.pyplot as plt
import geopandas as gpd
from shapely.geometry import LineString
from shapely.ops import split
from shapely.affinity import translate

def shift_geom(shift, gdataframe, plotQ=False, extent=[-60,120, -50,40]):
    # shift: expect positive values for good results
    # crs: crs="EPSG:4326" (not anything else); this handles world geometries
    # this code is adapted from answer found in SO
    # will be credited here: ???
    shift -= 180
    moved_geom = []
    splitted_geom = []
    border = LineString([(shift,90),(shift,-90)])

    for row in gdataframe["geometry"]:
        splitted_geom.append(split(row, border))
    for element in splitted_geom:
        items = list(element)
        for item in items:
            minx, miny, maxx, maxy = item.bounds
            if minx >= shift:
                moved_geom.append(translate(item, xoff=-180-shift))
            else:
                moved_geom.append(translate(item, xoff=180-shift))

    # got `moved_geom` as the moved geometry
    # must specify CRS here
    moved_geom_gdf = gpd.GeoDataFrame({"geometry": moved_geom}, crs="EPSG:4326")

    # can change crs here
    if plotQ:
        #fig1, ax1 = plt.subplots(figsize=[8,6])
        fig = plt.figure( figsize=(8,8) )
        ax1 = fig.add_subplot( projection=ccrs.PlateCarree() )
        moved_geom_gdf.plot( ax=ax1, figsize=(8,5), scheme='quantiles', cmap='tab20')
        ax1.set_extent(extent, crs=ccrs.PlateCarree())
        # ax1.gridlines( draw_labels=True ) # careful, misleading labels
        plt.show()

    return moved_geom_gdf

wideseas = geopandas.GeoDataFrame.from_file('./data/Longhurst_world_v4_2010.shp')
wideseas180 = shift_geom(180, wideseas, True, extent=[-60,120, -50,70])

输出图：

长呼情节

于 2020-12-03T03:58:00.247 回答
0

您可以使用以下命令获得基础形状POLYGONS和MULTIPOLYGON（有趣的是两者都有）的坐标。

provinces['geometry']

在 geopandas 中，几何图形始终存储在['geometry']列中。您还可以使用内置方法：

provinces.geometry

或对于特定行：

provinces.geometry[i]

还要看crs

provinces.crs

返回

<Geographic 2D CRS: EPSG:4326>
Name: WGS 84
Axis Info [ellipsoidal]:
- Lat[north]: Geodetic latitude (degree)
- Lon[east]: Geodetic longitude (degree)
Area of Use:
- name: World
- bounds: (-180.0, -90.0, 180.0, 90.0)
Datum: World Geodetic System 1984
- Ellipsoid: WGS 84
- Prime Meridian: Greenwich

也许您可以使用geopandas 的投影管理工具从这里重新投影。

于 2020-12-02T13:24:28.990 回答