检查点是否属于shapefile（多边形和点在不同的单元中）

发布于 2025-02-06 14:32:49 字数 3080 浏览 2 评论 0原文

我正在尝试创建一个口罩，遵循问题的答案，以便我只能在形状文件中为区域上色。我将在此处发布代码示例（在我的问题中的代码几乎相同，但其中包括解决方案代码）。

我的代码的快速摘要：这是（LON，LAT，TEM）的数据。我将LON和LAT转换为ccrs.platecaree（）坐标系统（该系统为米而不是度数），然后将转换的坐标放入XP，YP中。然后将3个变量XP，YP和TEM插值，然后使用PMESHCOLOR（）绘制颜色映射（我上一个问题中与上面链接的相同颜色映射）。现在，我要限制颜色映射仅在形状文件内部（这是灰色线）。代码如下：

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.colors as clr
import pandas as pd
import numpy as np

from metpy.interpolate import interpolate_to_grid, remove_nan_observations
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

canada_east = -95
canada_west = -101.8
canada_north = 52.8
canada_south = 48.85

central_lon = (canada_east + canada_west)/2
central_lat = (canada_north + canada_south)/2

crs = ccrs.LambertConformal(central_longitude = central_lon, central_latitude = central_lat)
lat = np.array([49.8134 50.904  50.698  49.095  49.436  49.9607 49.9601 49.356  50.116

49.402  52.3472 50.411  49.24   49.876  49.591  49.905  49.498  49.088
49.118  50.5947 49.3776 49.148  49.1631 51.358  49.826  50.4324 49.96
49.68   49.875  50.829  51.572])
lon = np.array([-100.3721  -97.273   -99.068   -97.528  -100.308   -98.9054  -98.6367
-99.248   -96.434  -100.93   -101.1099 -100.893  -100.055   -99.909
-97.518   -99.354   -98.03    -99.325   -99.054   -98.0035 -100.5387
-100.491   -97.1454 -100.361   -96.776   -99.4392  -97.7463  -97.984
-95.92    -98.111  -100.488])
tem = np.array([-8.45   -4.026  -5.993  -3.68   -7.35   -7.421  -6.477  -8.03   -3.834
-13.04   -4.057  -8.79   -6.619 -10.89   -4.465  -8.41   -4.861  -9.93
-7.125  -4.424 -11.95   -9.56   -3.86   -7.17   -4.193  -7.653  -4.883
-5.631  -3.004  -4.738  -8.81])

xp, yp, _ = crs.transform_points(ccrs.PlateCarree(), lon, lat ).T
xp, yp, tem = remove_nan_observations(xp, yp, tem)

alt_x, alt_y, data = interpolate_to_grid( xp, yp, tem, minimum_neighbors=2, search_radius=240000, interp_type = 'barnes', hres = 1000)

# Create the figure and grid for subplots
fig = plt.figure(figsize=(17, 12))

# Main ax
ax = plt.subplot(111, projection=crs)
ax.set_extent([canada_west, canada_east, canada_south, canada_north], ccrs.PlateCarree())


# Read shape file again using geopandas and create a mask

cat_gdf = geopandas.read_file('MB_AGregion_Perim_South.shp')
cat_gdf = cat_gdf.to_crs(epsg = 4326)
mask = shapely.vectorized.contains(cat_gdf.dissolve().geometry.item(), alt_x, alt_y)

 print(mask)

使

[[False False False ... False False False]
 [False False False ... False False False]
 [False False False ... False False False]
 ...
 [False False False ... False False False]
 [False False False ... False False False]
 [False False False ... False False False]]

Process finished with exit code 0

我的alt_x和alt_y为米，我的形状文件几何形状在程度上。我相信这就是掩模包含所有虚假值（这意味着掩盖一切）的原因。但是，我找不到解决此问题的解决方案。如果有人可以帮助我解决这个问题，我将非常感谢。

原文

I'm trying to create a mask following this question's answer, so that I can color areas inside my shape file only. I will post the code example here (which is pretty much the same with the code in my previous question, but with the solution code included).

A quick summary of my code: It's a data of (lon, lat, tem). I transform the lon and lat into ccrs.PlateCaree() coordinate system (which is in meter and not degree anymore), and put the transformed coordinates in xp, yp. The 3 variables xp, yp, and tem are then interpolated and then I use pmeshcolor() to draw a color map (the same color map in my previous question that I left the link above). Now I want to limit the color map to be inside the shape file only (which is the grey lines). Code is as below:

import cartopy.crs as ccrs
import cartopy.feature as cfeature
import matplotlib.pyplot as plt
import matplotlib.colors as clr
import pandas as pd
import numpy as np

from metpy.interpolate import interpolate_to_grid, remove_nan_observations
from cartopy.io.shapereader import Reader
from cartopy.feature import ShapelyFeature

canada_east = -95
canada_west = -101.8
canada_north = 52.8
canada_south = 48.85

central_lon = (canada_east + canada_west)/2
central_lat = (canada_north + canada_south)/2

crs = ccrs.LambertConformal(central_longitude = central_lon, central_latitude = central_lat)
lat = np.array([49.8134 50.904  50.698  49.095  49.436  49.9607 49.9601 49.356  50.116

49.402  52.3472 50.411  49.24   49.876  49.591  49.905  49.498  49.088
49.118  50.5947 49.3776 49.148  49.1631 51.358  49.826  50.4324 49.96
49.68   49.875  50.829  51.572])
lon = np.array([-100.3721  -97.273   -99.068   -97.528  -100.308   -98.9054  -98.6367
-99.248   -96.434  -100.93   -101.1099 -100.893  -100.055   -99.909
-97.518   -99.354   -98.03    -99.325   -99.054   -98.0035 -100.5387
-100.491   -97.1454 -100.361   -96.776   -99.4392  -97.7463  -97.984
-95.92    -98.111  -100.488])
tem = np.array([-8.45   -4.026  -5.993  -3.68   -7.35   -7.421  -6.477  -8.03   -3.834
-13.04   -4.057  -8.79   -6.619 -10.89   -4.465  -8.41   -4.861  -9.93
-7.125  -4.424 -11.95   -9.56   -3.86   -7.17   -4.193  -7.653  -4.883
-5.631  -3.004  -4.738  -8.81])

xp, yp, _ = crs.transform_points(ccrs.PlateCarree(), lon, lat ).T
xp, yp, tem = remove_nan_observations(xp, yp, tem)

alt_x, alt_y, data = interpolate_to_grid( xp, yp, tem, minimum_neighbors=2, search_radius=240000, interp_type = 'barnes', hres = 1000)

# Create the figure and grid for subplots
fig = plt.figure(figsize=(17, 12))

# Main ax
ax = plt.subplot(111, projection=crs)
ax.set_extent([canada_west, canada_east, canada_south, canada_north], ccrs.PlateCarree())


# Read shape file again using geopandas and create a mask

cat_gdf = geopandas.read_file('MB_AGregion_Perim_South.shp')
cat_gdf = cat_gdf.to_crs(epsg = 4326)
mask = shapely.vectorized.contains(cat_gdf.dissolve().geometry.item(), alt_x, alt_y)

 print(mask)

which gives

[[False False False ... False False False]
 [False False False ... False False False]
 [False False False ... False False False]
 ...
 [False False False ... False False False]
 [False False False ... False False False]
 [False False False ... False False False]]

Process finished with exit code 0

My alt_x and alt_y is in meter and my shape file geometries are in degree. I believe this is the reason why the mask contains all false value (which means masking everything). However, I'm unable to find a solution to this problem. I would appreciate a lot if someone can help me with this problem.

分享到QQ

分享到微博

如果你对这篇内容有疑问，欢迎到本站社区发帖提问参与讨论，获取更多帮助，或者扫码二维码加入 Web 技术交流群。

发布评论

需要登录才能够评论，你可以免费注册一个本站的账号。

池木 2025-02-13 14:32:49

这是另一个MRE。保持事物非常简单 - 我只是使用（x，y）值的bbox和向量，而不是您拥有的网格。但这足以说明您要处理的问题。

import cartopy.crs as ccrs
import shapely.geometry
import shapely.vectorized
import numpy as np, pandas as pd, geopandas as gpd, matplotlib.pyplot as plt

canada_east = -95
canada_west = -101.8
canada_north = 52.8
canada_south = 48.85

central_lon = (canada_east + canada_west)/2
central_lat = (canada_north + canada_south)/2

crs = ccrs.LambertConformal(
    central_longitude=central_lon,
    central_latitude=central_lat,
)

x = np.arange(-120, -90)
y = np.ones_like(x) * 50

bbox = shapely.geometry.Polygon([
    (canada_west, canada_south),
    (canada_east, canada_south),
    (canada_east, canada_north),
    (canada_west, canada_north),
 ])

gdf = gpd.GeoDataFrame(geometry=[bbox], crs='epsg:4326')

这是X，Y和GeodataFrame的样子。请注意，它们全部是学位，Techincally WGS84，aka 'epsg：4326'：

In [2]: print(f"x values: {x}\n\ny values: {y}\n\ngeodataframe:\n{gdf}")
x values: [-120 -119 -118 -117 -116 -115 -114 -113 -112 -111 -110 -109 -108 -107
 -106 -105 -104 -103 -102 -101 -100  -99  -98  -97  -96  -95  -94  -93
  -92  -91]

y values: [50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
 50 50 50 50 50 50]

geodataframe:
                                            geometry
0  POLYGON ((-101.80000 48.85000, -95.00000 48.85...

如果我们使用Shapely.Dector.Dectorized.Contains 。：

In [3]: shapely.vectorized.contains(gdf.dissolve().geometry.item(), x, y)
Out[3]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False,  True,  True,  True,  True,  True,  True, False, False,
       False, False, False])

crs platecarree（）是一个绘制的crs，可以保留LAT/LON值，因此您可以将其视为等同于WGS84（只有差异是PlateCarree具有定义的映射到像素值）。因此，以下内容将我们的LAT/LON值转换为对象crs的CRS，它是lambertConformal：

In [4]: xp, yp, _ = crs.transform_points(ccrs.PlateCarree(), x, y ).T

现在，我们的xp和YP值以米为单位，看起来与GeodataFrame中的值有很大不同：

In [5]: print(f"xp values: {xp}\n\nyp values: {yp}\n\ngeodataframe:\n{gdf}")
xp values: [-1555398.30438114 -1484657.61843652 -1413737.15236547 -1342645.49406745
 -1271391.25217193 -1199983.05499594 -1128429.54949926 -1056739.40023735
  -984921.28831211  -912983.91032072  -840935.97730252  -768786.21368415
  -696543.35622316  -624216.15294999  -551813.36210868  -479343.75109632
  -406816.09540141  -334239.17754116  -261621.78599806  -188972.71415562
  -116300.7592336    -43614.72122269    29076.59818104   101764.396642
   174439.87225098   247094.22459092   319718.65580267   392304.37165025
   464842.58258582   537324.50481399]

yp values: [ 92537.27851361  75810.36405389  59862.8937992   44696.79886026
  30313.91572949  16715.98605863   3904.65644788  -8118.52175353
 -19352.09263517 -29794.69590175 -39445.06703778 -48302.03746074
 -56364.53466258 -63631.58233956 -70102.30051051 -75775.90562338
 -80651.71065011 -84729.12516982 -88007.65544033 -90486.90445793
 -92166.57200545 -93046.45468863 -93126.44596073 -92406.53613545
 -90886.8123881  -88567.45874503 -85448.75606135 -81531.08198695
 -76814.91092072 -71300.81395314]

geodataframe:
                                            geometry
0  POLYGON ((-101.80000 48.85000, -95.00000 48.85...

此时，如果我们使用包含在GeodataFrame中找到点，那么很明显，为什么我们只能获得false值：

In [6]: shapely.vectorized.contains(gdf.dissolve().geometry.item(), xp, yp)
Out[6]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False])

要将GeodataFrame转换为我们的目标crs，我们可以使用to_crs：

In [7]: gdf_lambert = gdf.to_crs(crs)

现在，我们的值阵容阵容，并且包含了预期的工作：

In [8]: print(f"xp values: {xp}\n\nyp values: {yp}\n\ngdf_lambert:\n{gdf_lambert}")
xp values: [-1555398.30438114 -1484657.61843652 -1413737.15236547 -1342645.49406745
 -1271391.25217193 -1199983.05499594 -1128429.54949926 -1056739.40023735
  -984921.28831211  -912983.91032072  -840935.97730252  -768786.21368415
  -696543.35622316  -624216.15294999  -551813.36210868  -479343.75109632
  -406816.09540141  -334239.17754116  -261621.78599806  -188972.71415562
  -116300.7592336    -43614.72122269    29076.59818104   101764.396642
   174439.87225098   247094.22459092   319718.65580267   392304.37165025
   464842.58258582   537324.50481399]

yp values: [ 92537.27851361  75810.36405389  59862.8937992   44696.79886026
  30313.91572949  16715.98605863   3904.65644788  -8118.52175353
 -19352.09263517 -29794.69590175 -39445.06703778 -48302.03746074
 -56364.53466258 -63631.58233956 -70102.30051051 -75775.90562338
 -80651.71065011 -84729.12516982 -88007.65544033 -90486.90445793
 -92166.57200545 -93046.45468863 -93126.44596073 -92406.53613545
 -90886.8123881  -88567.45874503 -85448.75606135 -81531.08198695
 -76814.91092072 -71300.81395314]

gdf_lambert:
                                            geometry
0  POLYGON ((-251935.076 -217891.780, 251935.076 ...

In [9]: shapely.vectorized.contains(gdf_lambert.dissolve().geometry.item(), xp, yp)
Out[9]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False,  True,  True,  True,  True,  True,  True,  True, False,
       False, False, False])

图形表示

我们也可以看到绘制时非常清楚。这是LAT/LON空间中的点和shapefile，使用包含掩码，以蓝色为蓝色，如果外部为蓝色，则在内部绿色：


In [12]: fig, ax = plt.subplots()
    ...: gdf.plot(ax=ax, alpha=0.5)
    ...: mask = shapely.vectorized.contains(gdf.dissolve().geometry.item(), x, y)
    ...: ax.scatter(x[mask], y[mask], color='green')
    ...: ax.scatter(x[~mask], y[~mask], color='blue')
Out[12]: <matplotlib.collections.PathCollection at 0x18f5a3ac0>

如果我们混合和匹配预测，则没有一个点处于形状（当轴通过米缩放时，这是靠近（0，0）的微小的不可见规格）：

In [14]: fig, ax = plt.subplots()
    ...: gdf.plot(ax=ax, alpha=0.5)
    ...: mask = shapely.vectorized.contains(gdf.dissolve().geometry.item(), xp, yp)
    ...: ax.scatter(xp[mask], yp[mask], color='green')
    ...: ax.scatter(xp[~mask], yp[~mask], color='blue')
Out[14]: <matplotlib.collections.PathCollection at 0x1919e58d0>

当Shapefile和要点都在LambertConformal中时，我们又重新开始营业：

In [16]: fig, ax = plt.subplots()
    ...: gdf_lambert.plot(ax=ax, alpha=0.5)
    ...: mask = shapely.vectorized.contains(gdf_lambert.dissolve().geometry.item(), xp, yp)
    ...: ax.scatter(xp[mask], yp[mask], color='green')
    ...: ax.scatter(xp[~mask], yp[~mask], color='blue')
Out[16]: <matplotlib.collections.PathCollection at 0x191a88ee0>

“ https://i.sstatic.net/j9cmw.png” alt =“弯曲点，有一些正确显示为相交的shapefile”>

Here's another MRE. Keeping things very simple - I'm just using the bbox and vector of (x, y) values rather than the mesh you have. But it should be enough to illustrate the issue you're dealing with.

import cartopy.crs as ccrs
import shapely.geometry
import shapely.vectorized
import numpy as np, pandas as pd, geopandas as gpd, matplotlib.pyplot as plt

canada_east = -95
canada_west = -101.8
canada_north = 52.8
canada_south = 48.85

central_lon = (canada_east + canada_west)/2
central_lat = (canada_north + canada_south)/2

crs = ccrs.LambertConformal(
    central_longitude=central_lon,
    central_latitude=central_lat,
)

x = np.arange(-120, -90)
y = np.ones_like(x) * 50

bbox = shapely.geometry.Polygon([
    (canada_west, canada_south),
    (canada_east, canada_south),
    (canada_east, canada_north),
    (canada_west, canada_north),
 ])

gdf = gpd.GeoDataFrame(geometry=[bbox], crs='epsg:4326')

Here's what x, y, and the geodataframe look like. Note that they're all in degrees, techincally WGS84, aka 'epsg:4326':

In [2]: print(f"x values: {x}\n\ny values: {y}\n\ngeodataframe:\n{gdf}")
x values: [-120 -119 -118 -117 -116 -115 -114 -113 -112 -111 -110 -109 -108 -107
 -106 -105 -104 -103 -102 -101 -100  -99  -98  -97  -96  -95  -94  -93
  -92  -91]

y values: [50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50 50
 50 50 50 50 50 50]

geodataframe:
                                            geometry
0  POLYGON ((-101.80000 48.85000, -95.00000 48.85...

If we use shapely.vectorized.contains when both the points and the geodataframe are in degrees, we get the expected answer:

In [3]: shapely.vectorized.contains(gdf.dissolve().geometry.item(), x, y)
Out[3]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False,  True,  True,  True,  True,  True,  True, False, False,
       False, False, False])

The CRS PlateCarree() is a plotting CRS which preserves lat/lon values, so you can think of it as equivalent to WGS84 (only difference is PlateCarree has a defined mapping to pixel values). So the following translates our lat/lon values into the CRS of our object CRS, which is LambertConformal:

In [4]: xp, yp, _ = crs.transform_points(ccrs.PlateCarree(), x, y ).T

Now, our xp and yp values are in meters and look very different from the values in the geodataframe:

In [5]: print(f"xp values: {xp}\n\nyp values: {yp}\n\ngeodataframe:\n{gdf}")
xp values: [-1555398.30438114 -1484657.61843652 -1413737.15236547 -1342645.49406745
 -1271391.25217193 -1199983.05499594 -1128429.54949926 -1056739.40023735
  -984921.28831211  -912983.91032072  -840935.97730252  -768786.21368415
  -696543.35622316  -624216.15294999  -551813.36210868  -479343.75109632
  -406816.09540141  -334239.17754116  -261621.78599806  -188972.71415562
  -116300.7592336    -43614.72122269    29076.59818104   101764.396642
   174439.87225098   247094.22459092   319718.65580267   392304.37165025
   464842.58258582   537324.50481399]

yp values: [ 92537.27851361  75810.36405389  59862.8937992   44696.79886026
  30313.91572949  16715.98605863   3904.65644788  -8118.52175353
 -19352.09263517 -29794.69590175 -39445.06703778 -48302.03746074
 -56364.53466258 -63631.58233956 -70102.30051051 -75775.90562338
 -80651.71065011 -84729.12516982 -88007.65544033 -90486.90445793
 -92166.57200545 -93046.45468863 -93126.44596073 -92406.53613545
 -90886.8123881  -88567.45874503 -85448.75606135 -81531.08198695
 -76814.91092072 -71300.81395314]

geodataframe:
                                            geometry
0  POLYGON ((-101.80000 48.85000, -95.00000 48.85...

At this point, if we use contains to find points within the geodataframe, it's clear why we get only False values:

In [6]: shapely.vectorized.contains(gdf.dissolve().geometry.item(), xp, yp)
Out[6]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False, False, False])

To convert the geodataframe to our target crs, we can use to_crs:

In [7]: gdf_lambert = gdf.to_crs(crs)

Now, our values line up, and contains works as expected once again:

In [8]: print(f"xp values: {xp}\n\nyp values: {yp}\n\ngdf_lambert:\n{gdf_lambert}")
xp values: [-1555398.30438114 -1484657.61843652 -1413737.15236547 -1342645.49406745
 -1271391.25217193 -1199983.05499594 -1128429.54949926 -1056739.40023735
  -984921.28831211  -912983.91032072  -840935.97730252  -768786.21368415
  -696543.35622316  -624216.15294999  -551813.36210868  -479343.75109632
  -406816.09540141  -334239.17754116  -261621.78599806  -188972.71415562
  -116300.7592336    -43614.72122269    29076.59818104   101764.396642
   174439.87225098   247094.22459092   319718.65580267   392304.37165025
   464842.58258582   537324.50481399]

yp values: [ 92537.27851361  75810.36405389  59862.8937992   44696.79886026
  30313.91572949  16715.98605863   3904.65644788  -8118.52175353
 -19352.09263517 -29794.69590175 -39445.06703778 -48302.03746074
 -56364.53466258 -63631.58233956 -70102.30051051 -75775.90562338
 -80651.71065011 -84729.12516982 -88007.65544033 -90486.90445793
 -92166.57200545 -93046.45468863 -93126.44596073 -92406.53613545
 -90886.8123881  -88567.45874503 -85448.75606135 -81531.08198695
 -76814.91092072 -71300.81395314]

gdf_lambert:
                                            geometry
0  POLYGON ((-251935.076 -217891.780, 251935.076 ...

In [9]: shapely.vectorized.contains(gdf_lambert.dissolve().geometry.item(), xp, yp)
Out[9]:
array([False, False, False, False, False, False, False, False, False,
       False, False, False, False, False, False, False, False, False,
       False,  True,  True,  True,  True,  True,  True,  True, False,
       False, False, False])

Graphical representation

We can also see this very clearly when plotting. Here's the points and shapefile in lat/lon space, using the contains mask to color the points blue if outside and green if inside:


In [12]: fig, ax = plt.subplots()
    ...: gdf.plot(ax=ax, alpha=0.5)
    ...: mask = shapely.vectorized.contains(gdf.dissolve().geometry.item(), x, y)
    ...: ax.scatter(x[mask], y[mask], color='green')
    ...: ax.scatter(x[~mask], y[~mask], color='blue')
Out[12]: <matplotlib.collections.PathCollection at 0x18f5a3ac0>

If we mix and match projections, none of the points are in the shape (which is a tiny invisible spec near (0, 0) when the axes are scaled by meters):

In [14]: fig, ax = plt.subplots()
    ...: gdf.plot(ax=ax, alpha=0.5)
    ...: mask = shapely.vectorized.contains(gdf.dissolve().geometry.item(), xp, yp)
    ...: ax.scatter(xp[mask], yp[mask], color='green')
    ...: ax.scatter(xp[~mask], yp[~mask], color='blue')
Out[14]: <matplotlib.collections.PathCollection at 0x1919e58d0>

When the shapefile and the points are both in LambertConformal, we're back in business:

In [16]: fig, ax = plt.subplots()
    ...: gdf_lambert.plot(ax=ax, alpha=0.5)
    ...: mask = shapely.vectorized.contains(gdf_lambert.dissolve().geometry.item(), xp, yp)
    ...: ax.scatter(xp[mask], yp[mask], color='green')
    ...: ax.scatter(xp[~mask], yp[~mask], color='blue')
Out[16]: <matplotlib.collections.PathCollection at 0x191a88ee0>