没戳，从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 数据汇总

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

可用的投影方式有：

注：本人自己翻译，有的没见过，可能不准，图片来源于上述链接。

通常，全世界的地图都是墨卡托（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个圆的信息

我梦寐以求的小公鸡，快点出来吧，这是用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！！！