努力使用cartopy和matplotlib绘制数据

问题描述 投票：0回答：1

我正在努力使用 cartopy （和 matplotlib）在底图顶部绘制点。我认为问题是我无法设法获取数据和底图以匹配轴的投影/crs，因为我的数据结构与我发现的大多数示例不同。不幸的是，其中一些（我的代码的顶部）无法更改，因为这意味着重写许多其他代码和函数：/

首先我尝试使用 matplotlib 绘制点并使用 cartopy 添加底图（我的代码中的版本 1）。根据我为 

ax
设置投影的位置，底图根本不会绘制，或者在没有所需扩展的情况下绘制在所有内容之上（-> 显示整个地球）。如果我按照建议添加
transform=
参数here，则根本不会绘制任何内容（至少在视口内不会绘制）。

在 geopandas.org 上，我找到了一个使用 cartopy 函数 (

add_geometries
) 添加点的示例（版本 2）。这会引发错误
Source CRS must be an instance of CRS or one of its subclasses, or None.
，我不知道它来自哪里。

这是我使用的代码。不幸的是，它相当冗长，但我不知道在不丢失必要信息的情况下我还能省略什么。

import pandas as pd                            #2.1.4
import geopandas as gpd                        #0.12.2

import matplotlib.pyplot as plt                #3.8.0
from shapely.geometry import Point, Polygon    #1.8.5

import cartopy.crs as ccrs                     #0.22.0
import cartopy.io.img_tiles as cimgt

#------------------------------------------
#-----------Creating data-environment------
#------------------------------------------

#------Part 1: has to stay--------
gdf_test1 = gpd.GeoDataFrame(
    {
        "geometry": [Point(13.10157, 52.36640), Point(13.10130, 52.36656), Point(13.10254, 52.36582), Point(13.10309, 52.36549), Point(13.10332, 52.36535)]
    },
    crs = 'EPSG:4326'
)

gdf_test2 = gpd.GeoDataFrame(
    {
        "geometry": [Point(13.43562, 52.48038), Point(13.50074, 52.45651), Point(13.50061, 52.45657), Point(13.49998, 52.45688), Point(13.49984, 52.45695)]
    },
    crs = 'EPSG:4326'
)

outline_test1 = Polygon((Point(12.984381, 52.354036), Point(13.147663, 52.354036), Point(13.147663, 52.412748), Point(12.984381, 52.412748), Point(12.984381, 52.354036)))
outline_test2 = Polygon((Point(13.62631009242705, 52.517735730410585) , Point(13.625180130891495, 52.50452520302402) , Point(13.621934955922162, 52.49144577310617) , Point(13.616607639378886, 52.478623231164505) , Point(13.609251157772427, 52.46618078976451) , Point(13.599937818411506, 52.45423790766038) , Point(13.5887585076519, 52.44290915361389) , Point(13.57582177055312, 52.432303120286605) , Point(13.561252732467345, 52.42252139799845) , Point(13.545191874144964, 52.413657617485086) , Point(13.527793672849741, 52.40579657006711) , Point(13.509225122743361, 52.399013412873614) , Point(13.489664148436754, 52.39337296594993) , Point(13.469297926125792, 52.388929107231334) , Point(13.448321127145887, 52.38572427048734) , Point(13.426934099105035, 52.3837890504421) , Point(13.405341000000002, 52.38314191835826) , Point(13.383747900894965, 52.3837890504421) , Point(13.362360872854113, 52.38572427048734) , Point(13.341384073874208, 52.388929107231334) , Point(13.321017851563246, 52.39337296594993) , Point(13.301456877256639, 52.399013412873614) , Point(13.282888327150259, 52.40579657006711) , Point(13.265490125855036, 52.413657617485086) , Point(13.249429267532655, 52.42252139799845) , Point(13.23486022944688, 52.432303120286605) , Point(13.221923492348102, 52.44290915361389) , Point(13.210744181588494, 52.45423790766038) , Point(13.201430842227573, 52.46618078976451) , Point(13.194074360621114, 52.478623231164505) , Point(13.188747044077843, 52.49144577310617) , Point(13.185501869108505, 52.50452520302403) , Point(13.18437190757295, 52.517735730410585) , Point(13.185369938647879, 52.53095019147478) , Point(13.188488252991643, 52.5440412712636) , Point(13.193698653775746, 52.5568827315916) , Point(13.200952657395336, 52.569350632900345) , Point(13.210181894677728, 52.581324538061054) , Point(13.221298711298406, 52.59268868614871) , Point(13.234196963913654, 52.60333312435861) , Point(13.248753006259824, 52.61315478651281) , Point(13.264826857184044, 52.62205850701603) , Point(13.282263540302736, 52.62995795966922) , Point(13.300894582770653, 52.63677651143157) , Point(13.320539658531636, 52.6424479820343) , Point(13.341008359453616, 52.64691730128411) , Point(13.362102075970094, 52.65014105694213) , Point(13.383615967296576, 52.65208792721094) , Point(13.405341000000002, 52.652738993096214) , Point(13.427066032703426, 52.65208792721094) , Point(13.448579924029907, 52.65014105694215) , Point(13.469673640546384, 52.64691730128411) , Point(13.490142341468363, 52.642447982034305) , Point(13.509787417229353, 52.63677651143158) , Point(13.528418459697264, 52.62995795966922) , Point(13.54585514281595, 52.62205850701603) , Point(13.56192899374018, 52.61315478651281) , Point(13.57648503608635, 52.60333312435861) , Point(13.589383288701594, 52.592688686148705) , Point(13.600500105322267, 52.581324538061054) , Point(13.609729342604664, 52.56935063290036) , Point(13.616983346224249, 52.5568827315916) , Point(13.622193747008357, 52.5440412712636) , Point(13.625312061352123, 52.53095019147478) , Point(13.62631009242705, 52.517735730410585)))

liste_test = pd.DataFrame({
    'order': [1, 2],
    'bez': ['Region 1 (circ)', 'Region 2 (rect)'],
    'gdf': [gdf_test1, gdf_test2],
    'filter': [outline_test1, outline_test2]
})

#------Part 2: could be changed--------

gdf_lang_test = gpd.GeoDataFrame(columns=['bez', 'geometry'], geometry='geometry', crs='EPSG:4326')  
for i, row in liste_test.iterrows():
    for j, row2 in row['gdf'].iterrows():
        ptX = row2['geometry'].x
        ptY = row2['geometry'].y
        gdf_this_test = gpd.GeoDataFrame([(row['bez'], Point(ptX, ptY))], columns=['bez', 'geometry'], geometry='geometry', crs='EPSG:4326')        #jeder Messwert innerhalb eines gdf
        gdf_lang_test = gpd.GeoDataFrame(pd.concat([gdf_lang_test, gdf_this_test], ignore_index=True))  #wird an neuen gdf angehängt

bland_test1 = gpd.read_file('https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/2_bundeslaender/1_sehr_hoch.geo.json')
bland_test2 = gpd.read_file('https://raw.githubusercontent.com/isellsoap/deutschlandGeoJSON/main/4_kreise/1_sehr_hoch.geo.json')

args = [bland_test1, bland_test2, gdf_lang_test]

#-----------------------------------------------
#-------------Actual Code-----------------------
#-----------------------------------------------

fig, ax = plt.subplots(figsize=(15, 15))
bland_crs = args[0].crs   #Returns: epsg:4326  (lower case!?)
data_crs = args[2].crs    #Returns: EPSG:4326  (upper case!?)

ax = plt.axes(projection=ccrs.Mercator())  

#-----PLOTTING DATA-----------------------------
#-Version 1----------
args[0].plot(ax=ax,facecolor="none", edgecolor='black', linewidths=1, alpha=1, figsize=(15,15))      #, transform=ccrs.PlateCarree())
args[1].plot(ax=ax,facecolor="none", edgecolor='grey', linewidths=0.5, alpha=1, figsize=(15,15))     #, transform=ccrs.PlateCarree()(
args[2].plot(ax=ax, zorder=2, markersize=6)     #, transform=ccrs.PlateCarree())

#-Version 2----------
#ax.add_geometries(args[0]['geometry'], facecolor="none", edgecolor='black', crs=bland_crs)
#ax.add_geometries(args[0]['geometry'], facecolor="none", edgecolor='black', crs=bland_crs)
#ax.add_geometries(args[2]['geometry'], facecolor="none", edgecolor='black', crs=data_crs)

#-----PLOTTING Outlines-------------------------
n=0
for i, row in liste_test.iterrows():
    n=n+1
    polygon_temp = row['filter']
    p = gpd.GeoSeries(polygon_temp, crs=data_crs)      
    args += (p,)                             
    #-Version 1----
    args[n+2].plot(ax=ax, facecolor="none")    #,transform=ccrs.PlateCarree())

    #-Version 2---
    #ax.add_geometries(args[n+2], facecolor="none", edgecolor='black', crs=data_crs)

#----Calculate Extent---------------------------
xmin, ymin, xmax, ymax = args[2].total_bounds
xmin = xmin*0.998
xmax = xmax*1.002
ymin = ymin*0.998
ymax = ymax*1.002
xlim = ([xmin, xmax])
ylim = ([ymin, ymax])
ax.set_xlim(xlim)
ax.set_ylim(ylim)

#----PLOTTING Background image--------------
#ax = plt.axes(projection=ccrs.Mercator())     #if uncommented, the basemap is drawn on top of everything with max extent (whole earth)

url = 'https://server.arcgisonline.com/ArcGIS/rest/services/World_Shaded_Relief/MapServer/tile/{z}/{y}/{x}.jpg'
ax.set_extent(ax.get_extent(), ccrs.Mercator())
image = cimgt.GoogleTiles(url=url)
ax.add_image(image, 1)

plt.show()

非常感谢您的帮助！ :/

python pandas matplotlib cartopy
1个回答
0
投票

我找到了一种显示大部分数据的方法。这并不完全是我想要的，因为情节是在 PlateCarree 而不是 Mercator 中，但现在对我来说很好，所以我就这样吧。

我将绘图的投影设置移至其自身的初始化，并稍微更改了添加底图和数据的顺序。代码现在如下所示：

#-----------------------------------------------
#-------------Actual Code-----------------------
#-----------------------------------------------

fig, ax = plt.subplots(figsize=(15, 15), subplot_kw={'projection': ccrs.PlateCarree()})
bland_crs = args[0].crs   #Returns: epsg:4326
data_crs = args[2].crs    #Returns: EPSG:4326

#----Calculate Extent---------------------------
xmin, ymin, xmax, ymax = args[2].total_bounds
xmin = xmin*0.998
xmax = xmax*1.002
ymin = ymin*0.998
ymax = ymax*1.002

ax.set_extent([xmin, xmax, ymin, ymax], ccrs.PlateCarree())
#----PLOTTING Background image--------------
#ax = plt.axes(projection=ccrs.Mercator())     #if uncommented, the basemap is drawn on top of everything with max extent (whole earth)
url = 'https://server.arcgisonline.com/arcgis/rest/services/Canvas/World_Light_Gray_Base/MapServer/tile/{z}/{y}/{x}.jpg'

image = cimgt.GoogleTiles(url=url)
ax.add_image(image, 13, interpolation='hanning', regrid_shape=3000)

#-----PLOTTING DATA-----------------------------
#-Version 1----------
#args[0].plot(ax=ax,facecolor="none", edgecolor='black', linewidths=1, alpha=1, figsize=(15,15))      #, transform=ccrs.PlateCarree())
#args[1].plot(ax=ax,facecolor="none", edgecolor='grey', linewidths=0.5, alpha=1, figsize=(15,15))     #, transform=ccrs.PlateCarree()(
args[2].plot(ax=ax, zorder=2, markersize=6)     #, transform=ccrs.PlateCarree())


#-----PLOTTING Outlines-------------------------
n=0
for i, row in liste_test.iterrows():
    n=n+1
    polygon_temp = row['filter']
    p = gpd.GeoSeries(polygon_temp, crs=data_crs)      
    args += (p,)                             
    #-Version 1----
    args[n+2].plot(ax=ax, facecolor="none")    #,transform=ccrs.PlateCarree())


plt.show()

我还在 

interpolation='hanning', regrid_shape=3000
中添加了 
add_image
，以获得更好的图像质量，如此处推荐。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.