用matlab修改WRF静态土地利用数据

心动勿行动

楼主WRF菜鸟一枚，这两天在倒腾怎样修改土地利用数据，将过程发一下，请大神们批评指正

总的来说，修改USGS原始数据的landuse可分为三步：a 将所修改的新数据写成binary形式 b 为数据创建一个索引文件 c 修改GEOGRID.TBL，以下做具体说明。

使用matlab进行修改，该代码并不是我写的，是在网上找的，源代码网址http://hydrometeorology.princeton.edu/wiki/index.php/WRF_-_Modify_USGS_landuse，该代码针对30s的数据进行修改，而我要修改的是15s数据，在源代码的基础之上稍加修改就完成了

1 我的原始数据存放目录 /data/WRF_INPUT_DATA/geog/modislanduse08_15s

查看 index ，主要内容如下（并非全部信息，但是包含了修改代码所需要的信息）

type=categorical

category_min=1

category_max=20

projection=regular_ll

dx=0.004166667 （数据精度为15s）

dy=0.004166667

known_x=1.0     

known_y=1.0

known_lat=10.00208333

known_lon=70.00208333

wordsize=1

tile_x=2400  (存放的格点数为2400x2400)

tile_y=2400

tile_z=1

2 根据上述网址中的代码针对以上index的内容进行修改，第一部分程序如下

% 1) 定义所要修改的区域范围

clc;

clear;

Leflon = 106.613;   % Veryfy this

Rightlon = 106.762; % Veryfy this

Bolat = 35.5031;     % Veryfy this

Toplat = 35.5949;    % Veryfy this

% 2) 根据上述范围，找到自己所要修改的数据具体存在于哪个文件中

% 文件的存储结构可参见http://www.mmm.ucar.edu/wrf/users/tutorial/201201/WPS-advanced.ppt.pdf

rows(1)=10.00208333; % 即为index中的known_lat

cols(1)=70.00208333; % 和 known_lon

size=2400;   

res=1/3600*15; %15s分辨率

coloc(1)=1;

roloc(1)=1;

for i=2:36

   cols(i)=cols(i-1)+size*res;

   coloc(i)=coloc(i-1)+size;

end

for i=2:18

   rows(i)=rows(i-1)+size*res;

   roloc(i)=roloc(i-1)+size;

end

temp1=find(rows<Bolat);

temp2=find(rows>Toplat);

for i=1:(temp2(1)-temp1(end))

   lat_st(i)=roloc(temp1(end)+(i-1));

end

new_ll_lat=rows((temp1(end)));

clear temp1 temp2;

temp1=find(cols<Leflon);

temp2=find(cols>Rightlon);

for i=1:(temp2(1)-temp1(end))

   lon_st(i)=roloc(temp1(end)+(i-1));

end

new_ll_lon=cols((temp1(end)));

clear temp1 temp2;

display('These are the tiles that you need to download into your own directory:');

for i=1:length(lon_st);

   for j=1:length(lat_st);

   display(strcat(num2str(lon_st(i)),'-*****.',num2str(lat_st(j)),'-*****'));

   end

end

执行此部分，得到如下信息

These are the tiles that you need to download into your own directory:

7201-*****.4801-*****

3 从linux系统中下载该数据，放到e:\matlab之下

说明：我要修改的区域比较小，只需要修改一个文件，对于修改范围稍大的区域可能会出现多个文件，源程序中就是两个文件，需要修改多个文件的参照源程序

4 根据需要修改数据。我的目的是将所要修改区域内的LU=10（草地）改为LU=12（农田）

fid1 = fopen ('e:\matlab\07201-09600.04801-07200.backup','rb');

Bdata=fread(fid1)

fclose(fid1);

Adata=reshape(Bdata,2400,2400);

Adata=rot90(Adata); %旋转数据是由数据存储结构决定的

new_cols(1)=new_ll_lon;

new_rows(1)=new_ll_lat;

for i=2:2400

   new_cols(i)=new_cols(i-1)+res;

end

for i=2:2400

   new_rows(i)=new_rows(i-1)+res;

end

temp1=find(new_rows<Bolat);

temp2=find(new_rows>Toplat);

lat_st=temp1(end);

lat_en=temp2(1);

clear temp1 temp2;

temp1=find(new_cols<Leflon);

temp2=find(new_cols>Rightlon);

lon_st=temp1(end); % important

lon_en=temp2(1);   % important

clear temp1 temp2;

%修改数据

for i=lat_st:lat_en

   for j=lon_st:lon_en

       if Adata(i,j)==10 % 草地

           Adata(i,j)=12; % 农田

       end

   end

end

%可以将数据输出到一个txt文件中，检查是否为想要的数据，该步骤非必需

fid=fopen('e:\matlab\modify_test.txt','w');

fprintf(fid,'%s\n','ncols 2400');

fprintf(fid,'%s\n','nrows 2400');          % 根据 Adata修改

fprintf(fid,'%s\n','xllcorner 100.00208'); % 根据 new_ll_lon修改

fprintf(fid,'%s\n','yllcorner  30.00208'); % 根据new_ll_lat修改

fprintf(fid,'%s\n','cellsize 0.00416667');

fprintf(fid,'%s\n','NODATA_value -99.000000');   

for i=1:2400                                    

   for j=1:2400

         fprintf(fid,'%6.2f',Adata(i,j));

   end

         fprintf(fid,'\n');

end

fclose(fid);

%最终目的：将数据写成二进制格式

fid2=fopen('e:\matlab\07201-09600.04801-07200','wb');

topnyc=rot90(rot90(rot90(Adata(1:2400,:))));

count11 = fwrite(fid2,topnyc);

fclose(fid2);

5 到此程序修改完毕，得到一个07201-09600.04801-07200的二进制数据文件，接下来要为新数据创建索引。

在linux系统下建立一个新目录（我的为/home/macc/data/PL） 将/原来存放数据的目录（我的为data/WRF_INPUT_DATA/geog/modislanduse08_15s）下的所有文件（包括index和所有的数据文件）都拷到新目录下，然后将修改过后的新数据（在此是07201-09600.04801-07200）上传到/home/macc/data/PL，替换原来的07201-09600.04801-07200数据，因为所有的数据存储格式和数据数量都与原文件相同，因此index文件不需要修改

6 修改GEOGRID.TBL 中LANDUSEF部分的相应数据的路径，我的修改如下

将name=LANDUSEF部分

rel_path=     modis_15s08:modislanduse08_15s/  修改为

abs_path=      modis_15s08:/home/macc/data/PL

7 到此全部修改完毕，运行./geogrid.exe ，用ncview查看结果，与原数据对比，查看是否修改成功

由于楼主水平极为有限，错误在所难免，请看官批判着看，欢迎批评指正

pkmn99

  这一行代码有误： lon_st(i)=roloc(temp1(end)+(i-1)); roloc应该是coloc

pkmn99

我也写了一小段，功能是将多个区域拼起来，进而可以存成tif等其它格式
data=[];%save combined file
for i=1:length(lon_st);
    data=[data,temp];
    temp=[];
   for j=1:length(lat_st);
   name=strcat(num2str(lon_st(i)),'-',num2str(lon_st(i)+1199),'.',num2str(lat_st(j)), ...
   '-',num2str(lat_st(j)+1199)); %creat file name
   fid1 = fopen (['/User/Data/landuse_30s/',name],'rb');%data path
   Bdata=fread(fid1);
   fclose(fid1);
   Adata=reshape(Bdata,1200,1200);%1200, depend on size
   Adata=rot90(Adata);
   temp=[Adata;temp];  
   end
end

liutaoz21

顶楼主~~~~~~

kongfeng0824

感谢分享！！！！！！！！！

冰紫蓦然

顶一个

jinziqiqi1991

好东西  

climate05

论坛上很多人讨论如何修改静态数据，这篇帖子给出了方向和方法，应该顶上去

月贝凡

不错，支持            

147

收藏，之后使用方便

一颗永恒的流星

不错，学习了，{:eb502:}{:eb502:}