GeoHash在matlab中工具包编写

1.基础知识

以下是我之前看到不错的博文

1.GeoHash核心原理解析 http://www.cnblogs.com/LBSer/p/3310455.html
2.GeoHash距离大致估计 http://blog.csdn.net/bitcarmanlee/article/details/55824141
3.C实现GeoHash编码 http://blog.csdn.net/chen77716/article/details/5576668
4.基于GEOHASH算法的附近点搜索实现 http://www.cnblogs.com/subaiBlog/p/5410590.html

2.Matlab代码

1.对经纬度进行编码

function geohash=encode(lat,lot,prec) %编码% Base32编码    Base32= '0123456789bcdefghjkmnpqrstuvwxyz';latitude=[-90,90];%纬度longitude=[-180,180];%经度length=prec*5;%需要的二进制编码长度k=1;bits=zeros(1,5);%记录二进制码l=1;for i=1:length    if mod(i,2)~=0          mid=mean(longitude);        if lot>mid            bits(l)=1;            l=l+1;            longitude(1)=mid;        else            bits(l)=0;            l=l+1;            longitude(2)=mid;        end    else        mid=mean(latitude);  %从1开始奇数位为纬度        if lat>mid            bits(l)=1;            l=l+1;            latitude(1)=mid;        else            bits(l)=0;            l=l+1;            latitude(2)=mid;        end    end    if ~mod(i,5)        geohash(k)=Base32(1*bits(5)+2*bits(4)+4*bits(3)+8*bits(2)+16*bits(1)+1);        bits=zeros(1,5);        l=1;        k=k+1;    endend

2.将geohash编码解码成经纬度

function [lat,lot]=docode(geoh)      %解码% 编码长度    Len=size(geoh,2);% Base32编码    Base32= '0123456789bcdefghjkmnpqrstuvwxyz';% geoh为待解析的编码，是字符串%当前计算位的奇偶性odd = true ;latitude=[-90,90];%纬度longitude=[-180, 180];%经度for i=1:Len   for j=1:32            if(Base32(j)==geoh(i))                bits=j-1;    % 找到第i个字符对应的数                break;            end    end    for jj=1:5            j=6-jj;        switch j                case 5                ad=16;            case 4                ad=8;            case 3                ad=4;            case 2                ad=2;            case 1                ad=1;        end          if  bitand(bits,ad)              bit=1;          else              bit=0;          end      % 取出对应的位            if odd                mid = (longitude(1) + longitude(2)) / 2;                if bit==0                    longitude(2)= mid;                else                    longitude(1)= mid;                end            else                mid = (latitude(1) + latitude(2)) / 2;                if bit==0                    latitude(2)= mid;                else                    latitude(1)= mid;                end            end            odd = ~odd;    endend    lat = (latitude(1) + latitude(2)) / 2;    lot = (longitude(1) + longitude(2)) / 2;end

3.根据当前区域的GeoHash，推算出周围8个方位区域块的的GeoHash值。使用递归实现，输出为九个宫格形式。

function expand=geohash_expand(geoh)% 编码长度    Len=size(geoh,2);% Base32编码    Base32= '0123456789bcdefghjkmnpqrstuvwxyz';% geoh为待解析的编码，是字符串n=0;str1='';str='';expand=cell(3);for i=1:Len    n=find(Base32==geoh(i));    n=n-1;    str1=dec2bin(n,5);    str=[str,str1];endcode=str;len=size(code,2);              code=findup(code,len);            %九宫格，第一个geohashcode=findleft(code,len);expand(1,1)={bin2geohash(code,len)};code=str;code=findup(code,len);expand(1,2)={bin2geohash(code,len)}; %九宫格，第二个geohashcode=str;code=findup(code,len);code=findright(code,len);expand(1,3)={bin2geohash(code,len)};%九宫格，第三个geohashcode=str;code=findleft(code,len);expand(2,1)={bin2geohash(code,len)};%九宫格，第四个geohashexpand(2,2)={geoh};    %九宫格，第五个geohash，也是正中心位置code=str;code=findright(code,len);expand(2,3)={bin2geohash(code,len)};%九宫格，第六个geohashcode=str;code=finddown(code,len);code=findleft(code,len);expand(3,1)={bin2geohash(code,len)};%九宫格，第七个geohashcode=str;code=finddown(code,len);expand(3,2)={bin2geohash(code,len)};%九宫格，第八个geohashcode=str;code=finddown(code,len);code=findright(code,len);expand(3,3)={bin2geohash(code,len)};%九宫格，第九个geohash%%用于将二进制字符串转换成geohashfunction geohash=bin2geohash(str,len)Base32= '0123456789bcdefghjkmnpqrstuvwxyz';geohash='';str1='';len=(len/5);for i=1:len    str1=Base32(bin2dec(str((i-1)*5+1:i*5))+1);    geohash=[geohash,str1];endfunction code=findup(code,len)if len<=1    return;endlen=len-2;if code(len+2) == '1'    code=findup(code,len);    code(len+2)='0';else    code(len+2)='1';endfunction code=finddown(code,len)if len<=1    return;endif code(len)=='0'    len=len-2;    code=finddown(code,len);    code(len+2)='1';else    code(len)='0';endfunction code=findleft(code,len)if len<=1    return;endlen=len-1;if code(len) == '0'    len=len-1;    code=findleft(code,len);  %不要写成code=findright(code,len-1),浪费了一下午找原因    code(len+1) = '1';else    code(len) = '0';endfunction code=findright(code,len)if len<=1    return;endlen=len-1;if code(len) == '1'    len=len-1;    code=findright(code,len);     code(len+1) = '0';else    code(len) = '1';end

实例：

这里得到九宫格是和上图中成镜像关系，也就是一个是顺时针，一个为逆时针，其余均相同。

以上代码均放在GitHub托管，欢迎下载https://github.com/Wangxingzhi/geohash