UTM坐标的遥感图，Python终于把它扶正啦(手动狗头)

之前发布了两个大图，一个矢量不对位(这是个明显错误)，看红圈圈里的地方，很明显让我国的国土面积少了一大片（别问我遥感底图是啥，我不会告诉你的，毕竟成果还没公开）。

后来修正过，虽然没什么大问题，但总觉得怪怪的，emmmmmmmmm，错的地方都对上了，但整个图歪了！！！强迫症的我怎么能忍。

让我们来详细分析一下问题，原po的途径是利用basemap来画图，basemap的map函数中有许多参数，其中对图像影响较大的应该是投影的选择，它的默认投影选择项较少：

Value  Description
cea  Cylindrical Equal Area
mbtfpq  McBryde-Thomas Flat-Polar Quartic
aeqd  Azimuthal Equidistant
sinu  Sinusoidal
poly  Polyconic
omerc  Oblique Mercator
gnom  Gnomonic
moll  Mollweide
lcc  Lambert Conformal
tmerc  Transverse Mercator
nplaea  North-Polar Lambert Azimuthal
gall  Gall Stereographic Cylindrical
npaeqd  North-Polar Azimuthal Equidistant
mill  Miller Cylindrical
merc  Mercator
stere  Stereographic
eqdc  Equidistant Conic
rotpole  Rotated Pole
cyl  Cylindrical Equidistant
npstere  North-Polar Stereographic
spstere  South-Polar Stereographic
hammer  Hammer
geos  Geostationary
nsper  Near-Sided Perspective
eck4  Eckert IV
aea  Albers Equal Area
kav7  Kavrayskiy VII
spaeqd  South-Polar Azimuthal Equidistant
ortho  Orthographic
cass  Cassini-Soldner
vandg  van der Grinten
laea  Lambert Azimuthal Equal Area
splaea  South-Polar Lambert Azimuthal
robin  Robinson

其实正常情况下也够用，但我这样的确实是奇葩，属于不够用的类型。。。后来我发现其中有个参数：

rsphere radius of the sphere used to define map projection (default6370997 meters, close to the arithmetic mean radius of the earth). If given as a sequence, the first two elements are interpreted as the radii of the major and minor axes of an ellipsoid. Note: sometimes an ellipsoid is specified by the major axis and an inverse flattening parameter (if). The minor axis (b) can be computed from the major axis (a) and the inverse flattening parameter using the formula if = a/(a-b).

no_rot only used by oblique mercator. If settoTrue, the map projection coordinates will not be rotated totrue North. DefaultisFalse (projection coordinates are automatically rotated).

几经调整，还是无果。所以只能另辟蹊径。

需要说明的是，我不想再用basemap，所以我需要准备一个自用的矢量文件，用以替代basemap的底图矢量。先把结果输出吧：

聪明的你应该立刻就能发现不同了吧，这个图的显示方式是彻底利用UTM投影，包括范围的单位都改用米制。而且关键的是，终于不歪脖子了。

思路就是把遥感图转为二维数组，同时计算每个数组的位置坐标，并且叠加上矢量就好了。

我觉得关键是叠加矢量：

patches = [PolygonPatch(feature, edgecolor="red", facecolor="none", linewidth=2) for feature in features]
ax.add_collection(mpl.collections.PatchCollection(patches, match_original=True))

另一部分是计算坐标：

p1 = pyproj.Proj(init='epsg:32650')
for i in range(len(projlatlon)):
    coor[i,0],coor[i,1]= p1(projlatlon[i,0],projlatlon[i,1],inverse=True)

完整代码关注公众号，发送关键词 ‘UTM’

您受累。。。‘UTM’一定是大写字母

数据暂时不能分享，有类似问题的同学们可以参照着修改~

转发关注本公众号，会得到up主精神支持+1~ （搞笑）

微信号：一个有趣的灵魂

关注我们，了解更多