在我们日常生活中
尤其是农民
经常要计算土地的面积
那么
有什么简便的计算方法呢？
传统计算土地的亩数的公式是：1亩=60平方丈=6000平方尺=666.667平方米。
我们平常求亩数，一般用更简易的计算方法。公式是：长（米）×宽（米）×0.0015=亩数
如：长200米，宽150米的土地面积是：
200（米）×150（米）×0.0015=30000平方米×0.0015=45（亩）
长（米）×宽（米）÷667=亩数
如：一块地长60米，宽40米，计箅公式就是60X40÷667算出来的数就是亩数。
民间还有一个更实用的口决来计算：平方米换为亩，计算口诀为“加半左移三”。
如：1平方米＝0.0015亩，如128平方米等于多少亩？
计算方法是先用128加128的一半：128＋64＝192，再把小数点左移3位，即得出亩数为0.192。
用这些方法是不是计算起来要简便很多呢？赶紧收藏吧，以后遇到计算土地亩数的时候就可以用啦！

凸包模板题。。。
我用的是一种叫Graham的方法，主要思想是以最下最左的点作为原点，将其余各点按照极角排序（用向量的叉积实现），之后用栈维护新点永远在凸包最后两点的射线右边即可。
用反证法不难证明四个点一定全在凸包上（不然一定存在一个比已知四边形更大的四边形），枚举对角线将四边形隔成两个三角形，旋转卡壳维护两个三角形最大面积即可。
#include<cstdlib>
#include<cstdio>
#include<iostream>
#include<cstring>
#include<cmath>
#include<algorithm>
#include<vector>
#include<queue>
using namespace std;
#define maxn 2005
struct P{double x,y;}a[maxn],s[maxn];
int top,n;
double dis(P x,P y)
{
	return sqrt((x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y));
}
P operator +(P x,P y)
{
	return (P){x.x+y.x,x.y+y.y};
}
P operator -(P x,P y)
{
	return (P){x.x-y.x,x.y-y.y};
}
double operator *(P x,P y)
{
	return x.x*y.y-x.y*y.x;
}
bool operator < (P x,P y)
{
	double t=(x-a[1])*(y-a[1]);
	if(t==0)return dis(x,a[1])<dis(y,a[1]);
	return t>0;
}
void Init()
{
	scanf("%d",&n);
	for(int i=1;i<=n;i++)
	{
		scanf("%lf%lf",&a[i].x,&a[i].y);
	}
}
void graham()
{
	top=1;
	for(int i=2;i<=n;i++)if((a[i].x<a[top].x)||(a[i].x==a[top].x && a[i].y<a[top].y))top=i;
	if(top!=1)swap(a[1],a[top]);
	sort(a+2,a+1+n);
	top=2; s[1]=a[1]; s[2]=a[2];
	for(int i=3;i<=n;i++)
	{
		while(top>1 && (s[top]-s[top-1])*(a[i]-s[top-1])<=0)top--;
		s[++top]=a[i];
	}
	return ;
}
void work()
{
	double ans=0;
	int nc,nd;
	for(int i=1;i<=top;i++)
	{
		nc=i+1; nd=i+3;
		for(int j=i+2;j<=top;j++)
		{
			while(nc%top+1!=j && (s[nc%top+1]-s[nc])*(s[j]-s[i])>=0)
			{
				nc=nc%top+1;
			}
			while(nd%top+1 !=i && (s[nd%top+1]-s[nd])*(s[i]-s[j])>=0)
			{
				nd=nd%top+1;
			}
			ans=max(ans,(s[nc]-s[i])*(s[j]-s[i])+(s[j]-s[i])*(s[nd]-s[i]));
		}
	}
	ans=ans/2.0;
	printf("%.3lf\n",ans);
	return ;
}
int main()
{
	freopen("in.txt","r",stdin);
	Init();
	graham();
	work();
	return 0;
}

前言

假设有两期用地数据2012年和2017年，分成三种用地类型，首先得到这两期数据的栅格图像，每个栅格单元存储的值代表一种用地类型，然后以矩阵的方式读入内存中，再用numpy的flatten的方法将其一维化，最后调用python机器学习库sklearn.metrics中的confusion_metrix方法进行计算。

from sklearn.metrics import confusion_matrix
import numpy as np

a = np.array([[1, 2], [3, 1]])  # 代表2012用地数据
b = np.array([[1, 3], [1, 2]])  # 代表2017
a = a.flatten()
b = b.flatten()
res = confusion_matrix(a, b)
print(res)

#  输出结果：
#  [[1 1 0]
#  [0 0 1]
#  [1 0 0]]

取结果第一行第一列的值（1，1）= 1，表示2012年为1的用地类型在2017年也为1的用地类型有1个。（1， 2） = 1 表示2012年为1的用地类型在2017年也为2的有1个地块。

实际案例

以四川2000年和2010年土地利用变化数据为例，首先将两份数据保存为png格式的栅格图片，然后读入Python中，构建混淆矩阵，即可计算出土地利用面积转移矩阵：

import numpy as np
import sys, os
sys.path.append(os.pardir)
from PIL import Image
from sklearn.metrics import confusion_matrix


sc_lucc_2000_img = Image.open('F:\\Temp\\SCLUCC2000.png')
sc_lucc_2010_img = Image.open('F:\\Temp\\SCLUCC2010.png')
sc_lucc_2000_array = np.asarray(sc_lucc_2000_img).flatten()
sc_lucc_2010_array = np.asarray(sc_lucc_2010_img).flatten()
transition_area = confusion_matrix(sc_lucc_2000_array, sc_lucc_2010_array)
print(transition_area)

计算结果：



与另一博主（数据由他提供）计算结果完全一致：https://blog.csdn.net/XiaoQingwaShine/article/details/95501091?utm_medium=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromBaidu-1.control&depth_1-utm_source=distribute.pc_relevant_t0.none-task-blog-BlogCommendFromBaidu-1.control

但计算和操作过程简便了许多。

微信

欢迎关注我的微信鼓励我继续写文章分享知识~



 

 

1.叶面积指数（LAI）的定义：

叶面积指数，亦称叶面积系数。是指单位土地面积上植物叶片总面积占土地面积的倍数。它与植被的密度、结构（单层或复层）、树木的生物学特性（分枝角、叶着生角、耐荫性等）和环境条件（光照、水分、土壤营养状况）有关，是表示植被利用光能状况和冠层结构的一个综合指标。
2.经验模型预测LAI：

经验模型法基于数理统计方法，利用反射率、光谱一阶导数、植被指数等作为自变量，建立叶面积指数与特征光谱（如反射率、光谱指数）之间的回归模型，并利用该模型估算LAI。经验模型原理和算法比较简单，且估算精度满足要求，在植被LAI的估算中应用广泛，是最为常见的建模方式。
3.LAI经验模型原理（以三种下垫面为例，自己可根据实际修改）



4.程序源代码
;#############导入、输出文件及调用函数的程序 ############;
Pro LAI
;根据NDVI和Landcover计算LAI

;打开NDVI和Landcover文件
fn_NDVI=dialog_pickfile(title='选择NDVI数据',get_path=work_dir)
cd,work_dir
fn_landcover=dialog_pickfile(title='选择土地覆盖数据',get_path=work_dir)
cd,work_dir

;读入NDVI和Landcover数据
envi_open_file,fn_NDVI,r_fid = fid_NDVI
envi_open_file,fn_landcover,r_fid=fid_landcover
envi_file_query,fid_NDVI,ns=ns,nl=nl,nb=nb,dims=dims,$
data_type=data_type,interleave=interleave,offset=offset
map_info=envi_get_map_info(fid = fid_NDVI)
NDVI=envi_get_data(fid = fid_NDVI,dims=dims,pos=0)
landcover=envi_get_data(fid = fid_Landcover,dims=dims,pos=0)

;调用cal_LAI函数计算LAI
LAI_result=lai_cal(NDVI,landcover)

;保存结果
o_fn=dialog_pickfile(title='结果保存为')
envi_write_envi_file,LAI_result,out_name=o_fn,/no_copy,$
ns=ns,nl=nl,nb=nb,dims=dims,data_type=data_type,interleave=interleave,$
offset=offset,map_info=map_info
end

;#############叶面积指数计算程序############;
function lai_cal,ndvi,landcover
;计算 LAI 的函数
sz=size(ndvi)
result=make_array(size=sz)
;分段计算LAI
  w1=where(Landcover eq 1 and NDVI le 0.125,count)
  w2=where(Landcover eq 2 and NDVI ge 0.125 and NDVI le 0.825,count)
  w3=where(Landcover eq 2 and NDVI ge 0.825,count)
  w4=where(Landcover eq 3 and NDVI ge 0.125 and NDVI le 0.825,count)
  w5=where(Landcover eq 3 and NDVI ge 0.825,count)
if total(w1) gt 0 then result[w1]=0
if total(w2) gt 0 then result[w2]=0.1836*exp(4.37*NDVI[w2])
if total(w3) gt 0 then result[w3]=6.606
if total(w4) gt 0 then result[w4]=0.0884*exp(4.96*NDVI[w4])
if total(w5) gt 0 then result[w3]=6.091
return,result ;输出结果
end


5.调用过程：

1）在IDL中调用 ENVI 函数读取 NDVI 和土地覆盖图

2） cal_LAI函数按照不同土地覆盖类型的 LAI 计算公式进行计算

3） 将计算得到的 LAI 保存为 ENVI 标准文件
6.结果：

ENVI标准文件做彩色变换后，结果如下：