之前发布了两个大图,一个矢量不对位(这是个明显错误),看红圈圈里的地方,很明显让我国的国土面积少了一大片(别问我遥感底图是啥,我不会告诉你的,毕竟成果还没公开)。
后来修正过,虽然没什么大问题,但总觉得怪怪的,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~ (搞笑)
微信号:一个有趣的灵魂
关注我们,了解更多