碼上快樂

相關內容    簡體    繁體

Matlab 畫地圖之 m_map

本文轉載自   查看原文   2020-10-15 19:35   1772    Matlab

沒戳，從CSDN上搞過來的，今天（2020.10.15）一看，圖被和諧了，嚴重影響觀看效果。所以搞到這里來了。

參考：
氣象家園： https://mp.weixin.qq.com/s/JWi7VllWr93wDw6BStk8Dw
M-Map：https://www.eoas.ubc.ca/~rich/map.html
找shp的網站：https://gadm.org/download_country_v3.html
中國科學院資源環境數據平台：http://www.resdc.cn/Default.aspx
matlab利用m_map工具包畫中國地圖及散點雲圖：https://www.cnblogs.com/righdflf/p/11484189.html
百度經驗出圖：https://jingyan.baidu.com/article/870c6fc36fdacfb03ee4be58.html

快應用：
shp文件下載：
世界國家+中國大陸部分：https://download.csdn.net/download/Gou_Hailong/12744519
國家基礎地理數據：https://download.csdn.net/download/tessykandy/9112343
代碼下載：
求內切圓外接圓：https://download.csdn.net/download/Gou_Hailong/12744427
縮放矩陣或圖像：https://download.csdn.net/download/Gou_Hailong/12744441
扣shp：https://download.csdn.net/download/Gou_Hailong/12744543

@

目錄

1、shp數據

目前搞到的shape文件有：中國科學院資源環境數據、從CSDN上下載的、奇哥給的。

1.1 中國科學院資源環境數據

在這里插入圖片描述
這個數據精度很高，比如，數據2的P2 長度為925，這個數據P2長度是4461；當然精度高同時畫起來也超級慢，我的lj 電腦畫個圖就可以吃一碗泡面了。
它配套一個xls文件：中國各省份與地級市的編碼.xls里面存了省市代碼。P2數據存在下面的問題：有的省名字缺失，P3數據有的也缺失，但是其對應的代碼查詢不到，所以不管了。下面補一下P2的省名字。
在這里插入圖片描述
修正后的shp文件名為：sheng.shp。

1.2 從CSDN上下載的

下載鏈接：https://download.csdn.net/download/tessykandy/9112343
在這里插入圖片描述

L1file='D:\下載\Useful\shp\國家基礎地理數據\bou1_4m\bou1_4l.shp';%國界線
P1file='D:\下載\Useful\shp\國家基礎地理數據\bou1_4m\bou1_4p.shp';%國界多邊形
L2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4l.shp';%省界線
P2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4p.shp';%省界多邊形
P3file='D:\下載\Useful\shp\國家基礎地理數據\bou3_4m\diquJie_polyline.shp';%地市級
L4file='D:\下載\Useful\shp\國家基礎地理數據\bou4_4m\BOUNT_line.shp';%縣界線
P4file='D:\下載\Useful\shp\國家基礎地理數據\bou4_4m\BOUNT_poly.shp';%縣界多邊形
%讀取國界shp文件的內容
L1=shaperead(L1file);
P1=shaperead(P1file);
L2=shaperead(L2file);
P2=shaperead(P2file);
P3=shaperead(P3file);
L4=shaperead(L4file);
P4=shaperead(P4file);
bou1_4lx=[L1(:).X];%提取經度信息
bou1_4ly=[L1(:).Y];%提取緯度信息
bou1_4px=[P1(:).X];%提取經度信息
bou1_4py=[P1(:).Y];%提取緯度信息
bou2_4lx=[L2(:).X];%提取經度信息
bou2_4ly=[L2(:).Y];%提取緯度信息
bou2_4px=[P2(:).X];%提取經度信息
bou2_4py=[P2(:).Y];%提取緯度信息
bou3_4px=[P3(:).X];%提取經度信息
bou3_4py=[P3(:).Y];%提取緯度信息
bou4_4lx=[L4(:).X];%提取經度信息
bou4_4ly=[L4(:).Y];%提取緯度信息
bou4_4px=[P4(:).X];%提取經度信息
bou4_4py=[P4(:).Y];%提取緯度信息

figure;
subplot(231);
mapshow(L1,'Color','k');    title('1China mapshow PolyLine-1')
subplot(232);
mapshow(P1,'FaceColor','r');title('1China mapshow Polygon-2')
subplot(233);
geoshow(P2,'FaceColor','y');title('2China geoshow Polygon-3')
subplot(234);
mapshow(L2,'Color','g');    title('2China mapshow PolyLine-4')
subplot(235);
m_proj('miller','lon',[70,140],'lat',[10,60]);%選擇投影方式
m_plot(bou3_4px,bou3_4py);%繪國界, 默認黑色,
m_grid;%繪制邊框
title('3China m-proj PolyLine-5')
subplot(236);
m_proj('robinson','lon',[70,140],'lat',[10,60]);%選擇投影方式
m_plot(bou4_4lx,bou4_4ly,'b');%繪國界，藍色
m_grid;%繪制邊框
title('4China m-proj PolyLine-6')

結果：
在這里插入圖片描述
可以看到地級市只有內部線，想要使用，得和省界線配合使用

L2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4l.shp';%省界線
L3file='D:\下載\Useful\shp\國家基礎地理數據\bou3_4m\diquJie_polyline.shp';%地市級
L2=shaperead(L2file);
L3=shaperead(L3file);
bou2_4lx=[L2(:).X];%提取經度信息
bou2_4ly=[L2(:).Y];%提取緯度信息
bou3_4lx=[L3(:).X];%提取經度信息
bou3_4ly=[L3(:).Y];%提取緯度信息
m_plot(bou3_4lx,bou3_4ly,'k');
hold on;
m_plot(bou2_4lx,bou2_4ly,'linewidth',3,'color','k');

在這里插入圖片描述

1.3 奇哥給的

奇哥給的和上面那個是一樣的，多了世界國家。下面將兩個整合到一起。
在這里插入圖片描述
代碼用的是 3 中的功能1，然后將 NAME 改成 FCNAME 即可。調用：

Sjfile='D:\下載\Useful\shp\國家基礎地理數據\世界國家\世界國家.shp';
sj=shaperead(Sjfile);
drawsheng(Sjfile,'中國');

1.4 數據匯總

文件名 描述
bou1_4m 中國國界
bou2_4m 中國省界
bou3_4m 中國地市級（只有內部線）
bou4_4m 中國縣級
hyd1_4m 中國一級河流
hyd2_4m 中國二級河流
hyd4_4m 中國四級河流
hyd5_4m 中國五級河流
rai_4m 中國主要鐵路
res1_4m 省會標示
res2_4m 地市級以上居民地
res4_4m 縣以上居民地級
roa_4m 中國主要公路
世界國家 世界國家
省界 另一個省界（不用）
中國省級行政區划_shp.rar 一般不用
縣級_截止08年_shp.rar 一般不用
sheng 中國省界（精度高，慢）含修正后的
shi 中國市界（精度高，慢）

常用數據路徑匯總：

L1file='D:\下載\Useful\shp\國家基礎地理數據\bou1_4m\bou1_4l.shp';%國界線
P1file='D:\下載\Useful\shp\國家基礎地理數據\bou1_4m\bou1_4p.shp';%國界多邊形
L2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4l.shp';%省界線
P2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4p.shp';%省界多邊形
L3file='D:\下載\Useful\shp\國家基礎地理數據\bou3_4m\diquJie_polyline.shp';%地市級
L4file='D:\下載\Useful\shp\國家基礎地理數據\bou4_4m\BOUNT_line.shp';%縣界線
P4file='D:\下載\Useful\shp\國家基礎地理數據\bou4_4m\BOUNT_poly.shp';%縣界多邊形
P03file='D:\下載\Useful\shp\中國科學院資源環境數據\shi\CN-shi-A.shp';%地級市
P021file='D:\下載\Useful\shp\中國科學院資源環境數據\sheng\CN-sheng-A.shp';%省
P02file='D:\下載\Useful\shp\中國科學院資源環境數據\sheng\sheng.shp';%修正之后的省
Sjfile='D:\下載\Useful\shp\國家基礎地理數據\世界國家\世界國家.shp';

%讀取國界shp文件的內容
L1=shaperead(L1file);
P1=shaperead(P1file);
L2=shaperead(L2file);
P2=shaperead(P2file);
P02=shaperead(P02file);
P021=shaperead(P021file);
L3=shaperead(L3file);
P03=shaperead(P03file);
L4=shaperead(L4file);
P4=shaperead(P4file);
bou1_4lx=[L1(:).X];%提取經度信息
bou1_4ly=[L1(:).Y];%提取緯度信息
bou1_4px=[P1(:).X];%提取經度信息
bou1_4py=[P1(:).Y];%提取緯度信息
bou2_4lx=[L2(:).X];%提取經度信息
bou2_4ly=[L2(:).Y];%提取緯度信息
bou2_4px=[P2(:).X];%提取經度信息
bou2_4py=[P2(:).Y];%提取緯度信息
bou3_4lx=[L3(:).X];%提取經度信息
bou3_4ly=[L3(:).Y];%提取緯度信息
bou4_4lx=[L4(:).X];%提取經度信息
bou4_4ly=[L4(:).Y];%提取緯度信息
bou4_4px=[P4(:).X];%提取經度信息
bou4_4py=[P4(:).Y];%提取緯度信息

下面對中國行政單位划分做個總結：

https://www.cnblogs.com/Gou-Hailong/p/13545885.html
詳細描述請移駕：https://jingyan.baidu.com/article/59a015e398ed3ff794886589.html
行政區划代碼查詢：https://xingzhengquhua.51240.com/

2、m_map 使用筆記

2.1 安裝

下載鏈接：https://www.eoas.ubc.ca/~rich/map.html
在這里插入圖片描述

下好工具箱解壓后，將文件夾放在一個不常動的地方，我放在了D:\MATLAB\R2017a\toolbox中，然后打開matlab，點擊set path(設置路勁)，將m_map路徑加進去就行了。測試是否成功的話，就在matlab命令空間敲m_demo(1)括號里面可以是1-15，這實際上m_map中附帶的一個函數，在文件m_demo.m，看看可不可以畫出來圖。

2.2 常用命令

User Guide：https://www.eoas.ubc.ca/~rich/mapug.html#p2

m_proj get //得到當前地圖的投影方式
m_proj('set','xxx'); //xxx見下表，查看投影方式的細節
m_proj('xxx');       //設置當前投影方式
m_scale(250000);     //比例尺設成1:250000
m_scale('auto');     //默認比例尺
a=m_scale            //返回當前比例尺
m_coord('IGRF-geomagnetic',datenum(2000,1,1)); //設置IGRF-2000
m_coast('line', ...optional line arguments );
/*m_coast函數可以訪問M_Map數據庫。 
將海岸線繪制為簡單的線條。
可選參數是線的樣式、寬度、顏色等。
*/
m_coast('patch', ...optional patch arguments );
m_coast('patch',[.7 .7 .7],'edgecolor','g');//繪制灰色土地，綠色勾畫
m_coast('speckle', ....optional m_hatch arguments);

m_elev('contour',LEVELS,'edgecolor','b');//畫輪廓
m_elev('contourf',LEVELS, optional contourf arguments);//填充輪廓
[Z,LONG,LAT]=m_elev([LONG_MIN LONG_MAX LAT_MIN LAT_MAX]);//返回深度Z的矩形矩陣的位置,經緯度

m_plot(LONG,LAT,...line properties)      % draw a line on a map (erase current plot) 
m_plot(bou2_4lx,bou2_4ly,'linewidth',3,'color','k');
m_line(LONG,LAT,...line properties)      % draw a line on a map 
m_text(LONG,LAT,'string')                % Text   
m_quiver(LONG,LAT,U,V,S)                 % A quiver plot 
m_patch(LONG,LAT,..patch properties)     % Patches.  
m_annotation('line',LON,LAT)             % Annotation

可用的投影方式有：

參數 中文描述 圖
Stereographic 保角，常用於極區
Orthographic 既不是等面積也不是保角的，而是類似於地球的透視圖 在這里插入圖片描述
Azimuthal Equal-area 朗伯方位等面積投影 在這里插入圖片描述
Azimuthal Equidistant 等距
Gnomonic 點切投影，地圖上所有的直線(不僅僅是那些通過中心的)都是大圓路線
Satellite 地球的透視圖，由衛星在特定的高度所看到。指定視點高度而不是為地圖指定半徑
Albers Equal-Area Conic 等積阿波斯圓錐
Lambert Conformal Conic 蘭伯特等角圓錐 在這里插入圖片描述
Mercator 墨卡托，正形切圓柱
Miller Cylindrical 米勒圓柱 在這里插入圖片描述
Equidistant Cylindrical 等距圓柱
Cylindrical Equal-Area 等積圓柱
Oblique Mercator 斜軸墨卡托 在這里插入圖片描述
Transverse Mercator 橫軸墨卡托
Sinusoidal 偽圓柱投影
Gall-Peters 緯線和子午線的平行線都呈現為直線，但垂直尺度被扭曲
Hammer-Aitoff 具有彎曲子午線和平行線的等面積投影 在這里插入圖片描述
Mollweide 外德投影，也稱為橢圓或等積投影。在這個投影中，平行線是直的。 在這里插入圖片描述
Robinson 洛濱孫，養眼 在這里插入圖片描述
UTM UTM投影 在這里插入圖片描述

注：本人自己翻譯，有的沒見過，可能不准，圖片來源於上述鏈接。

通常，全世界的地圖都是墨卡托（Mercator），盡管米勒圓柱投影通常看起來更好，因為它沒有強調極地區域。 另一個選擇是Hammer-Aitoff或Mollweide（使子午線在兩極附近彎曲在一起）。 兩者都是等面積的。 將這些投影用於在中點附近沒有赤道的地圖可能不是一個好主意。 Robinson投影不是等面積或等角投影，而是國家地理雜志的選擇（無論如何），並且也出現在IPCC報告中。
如果您要繪制的東西的北/南范圍較大，但不是很寬（例如北美洲和南美洲，或者北大西洋和南大西洋），那么正弦或Mollweide投影將看起來非常不錯。 另一個選擇是“橫軸墨卡托”，盡管通常只用於非常大比例尺的地圖（即，覆蓋很小面積的地圖）。
對於一個半球或其他半球內較小的區域（例如，澳大利亞，美國，地中海，北大西洋），您可能會選擇圓錐投影。 兩種可用的圓錐投影之間的差異非常細微，如果您對投影不太了解，那么使用哪種圓錐投影可能不會有太大的區別。
如果您要畫的區域變得更小，則使用哪種投影都無關緊要。 我認為在很多情況下有用的一種投影是斜墨卡托，因為您可以將其沿着長（但窄）的沿海地區對齊。 如果沿着經度/緯度線的地圖限制正常，請使用橫軸墨卡托或圓錐投影。 如果要以米為單位建立網格，則UTM投影也很有用，因為投影坐標（即“東和北”）以米為單位。
極地地區傳統上是使用Stereographic投影來繪制地圖的，因為每個人似乎都認為具有“靶心”式的緯線模式看起來不錯（經度線在極點匯聚在一起的事實也使它們成為極地科學家幾乎不可抗拒的目的地 和探險者）。

3、常用功能實現

3.1 畫出某省並標出其在中國地圖上的位置

全代碼：https://blog.csdn.net/Gou_Hailong/article/details/108209764

調用實例：

file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4p.shp';%省界多邊形
str='廣東省';
drawsheng(file,str);
shapewrite(sheng,filename); //存到shp文件中
> 函數全代碼請移駕：[https://blog.csdn.net/Gou_Hailong/article/details/108209395](https://blog.csdn.net/Gou_Hailong/article/details/108209395)
![在這里插入圖片描述](https://img-blog.csdnimg.cn/20200821231423621.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0dvdV9IYWlsb25n,size_16,color_FFFFFF,t_70#pic_center)

3.2 用shp文件裁剪圖片

先搞個簡單的：河南省的
1.第一步，把河南省的shp提取出來，注意看其范圍

P2file='D:\下載\Useful\shp\國家基礎地理數據\bou2_4m\bou2_4p.shp';%省界多邊形
henan=drawsheng(P2file,'河南省');%第一步，把河南省的shp提取出來，注意看其范圍

可以看到，其大致范圍是110°E-117°E, 30°N-37°N。

2.從 tif 里面將該大致范圍搞出來Henan，並提取shp邊界xv, yv，之后用inpolygon搞出邏輯矩陣in(在邊界內的為1)

TDMSPfile='D:\復習資料\自學\7_科研\夜光遙感\data\DMSP\F121999.v4\F121999.v4b_web.stable_lights.avg_vis.tif';
TDMSP=imread(TDMSPfile);
Henan=TDMSP(38*120:45*120,290*120:297*120);%根據自己需要裁減
X=30:0.0083333333:37;
Y=110:0.0083333333:117;
[a,b]=size(Henan);
X1=repmat(X',1,b);
Y1=repmat(Y,a,1);
X1=flipud(X1);%這個得上下翻轉一下。
xv=henan.X;yv=henan.Y;        %提取邊界
xv = xv(1:end-1); yv = yv(1:end-1);
in = inpolygon(Y1,X1,xv,yv);    %返回邏輯矩陣，

3.把邊界外的數據搞成nan

photo=zeros(a,b);
photo(:)=nan;
for i=1:a
    for j=1:b
        if in(i,j)
            photo(i,j)=Henan(i,j);
        end
    end
end  %在邊界外的搞成nan

4.畫圖

[x,x1,y,y1] = getxy(Y,X);
x=(x-110).*120; 
y=(y-30).*120;
photo=flipud(photo);%contourf 上下翻轉一下，才變成imshow
contourf(photo,'LineStyle','none');
colormap(jet);colorbar
set(gca,'XTick',x,'XTicklabel',x1);   %設置x,y軸
set(gca,'YTick',y,'YTicklabel',y1);
title('河南省夜光遙感影像');

colormap(gray);colorbar  //寧夏的,

在這里插入圖片描述

PS：出邏輯矩陣的時候超級慢，這還是一個簡單的邊界，目前只能這么搞了，還沒有其他方法。
跑全中國的跑十幾分鍾還跑步不來。我丟，中國大陸矩陣大小4081x7561 = 30,856,441，本來以為是三百多萬，從8.23 14:50開始跑，跑到8.23 23:02跑了177萬，我覺得明天早上應該就跑完了，有點期待夾雜着莫名的興奮。到了8.24 7:57跑到了270萬，我想着在有倆小時估計就跑完了，但是10點之后，我發現還沒跑完，跑到了310萬。我去，不是308萬嗎，又仔細數了數零，發現是3080萬，當場哭暈，那不得10天嗎，就沒管它，它自己在那跑，到8.24 16:03 跑到390萬停，結果如下圖：
在這里插入圖片描述
算是把雄雞雞冠畫出來了吧。。

突然想到一個新方法，中國邊界內的最大內接圓一定在中國里面，中國邊界的最小外接圓一定在中國外邊，所以最后只需判斷中間夾的那個小圓環就好了。
在這里插入圖片描述
這是裁北京市的圖，小圓太小，大圓太大，效果還是不太理想，求重心求得精確一點就好了。

在這里插入圖片描述
改進了求重心的算法，效果如上圖（右）還可以。
這個功能的算法，目前來說，我是沒招了，全代碼鏈接如下：

https://blog.csdn.net/Gou_Hailong/article/details/108209764

但是，對於中國大陸來說，兩個圓還是有點少：
在這里插入圖片描述
上圖雖簡單，但上面花了我好幾個小時，搞出來了這7個圓的信息

R X Y 信息
11.807875718188898 1.338997219786234e+02 30.905277734156662 內部不要
31.334089276544790 1.044958332073500e+02 43.506944350416660 外部不要
13.553968551057642 1.338997219786233e+02 22.504166656650000 內部不要
11.586902594103240 1.044958332073500e+02 53.999999864000000 內部不要
9.333210635315217 1.086963887461033e+02 33.005555503533330 內部要
4.040350856671964 1.254986109011167e+02 45.607222119793335 內部要
8.597190190877908 87.693611052336680 37.206111042286665 內部要

我夢寐以求的小公雞，快點出來吧，這是用mapshow(shp)畫出來的。2020/8/24 20:34 還在苦逼等待程序運行中....（真累，明天換個活，先歇歇。）
2020/8/26 16:13 經歷了大約兩天時間，終於跑完了，加上那七個圓，算法運行時間大大減少了，之前的話少說也得10天，現在大約兩天跑好了，電腦燙的厲害。下面看下結果：

Elapsed time is 169864.909166 seconds.
169864.909166s = 47.18 h
------------1------------------
8.23 14.50 開始跑
8.23 20.05 1334000
8.23 23:02 1777000
8.24 7.57 2722000
8.24 14.30 374萬
8.24 16.03 390萬停
------------2------------------
8.24 14.22 開始跑
8.24 14.25 700多萬，慢起來了
8.24 14.30 766萬
8.24 16.02 1060萬
蹦了，程序錯了
------------3------------------
8.24 17:00 start
8.24 22:18 307w
8.25 6:45 713w
8.25 12:00 1000w
8.25 16:23 1239w, 做個備份，photo1  i=1638, j=6164
8.25 22:36 1905w, 做個備份，photo2  i=2520, j=664
8.26 8:53 2737w, 做個備份，photo3  i=3620, j=3749
8.26 16:03 ok; 
--------歷時-47.18 h--------

在這里插入圖片描述
這張圖片具有重大的“歷史紀念”意義，希望CSDN不要和諧掉。。
Amazing！！！


免責聲明！

本站轉載的文章為個人學習借鑒使用，本站對版權不負任何法律責任。如果侵犯了您的隱私權益，請聯系本站郵箱yoyou2525@163.com刪除。



猜您在找 Matlab——m_map指南（2） MATLAB——m_map指南（1） matlab利用m_map工具包畫中國地圖及散點雲圖 M_map(五) MATLAB M_map新增的函數可以添加指北針啦！ Matlab 畫地圖時搞定經緯度注釋 函數 matlab畫甘特圖 MATLAB之m文件的使用（1） matlab畫柱狀圖 用matlab畫漂亮的sin曲線
 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM