Showing posts with label ANUSPLIN. Show all posts

Sunday, October 15, 2023

Batch interpolation of meteorological variables with ANUSPLIN

ANUSPLIN is a globally renowned software for interpolating meteorological variables and has been widely used in meteorological interpolation activities around the world. In this instance, it is used to interpolate meteorological station data to obtain the temperature distribution for a particular area from December 1964 to 2022. The interpolation process utilizes Matlab code to coordinate and manage ANUSPLIN's input, operation, and output phases.

References

[1] ANUSPLIN.

Saturday, August 12, 2023

Manitoba: Analysis for the Month of February

As can be seen from the following figures, the February temperature in Manitoba has generally maintained a pattern of being higher in the south and lower in the north over 109 years. Upon closer examination of February 2023, unusually low temperatures were observed in the northeastern region, a situation not seen in many years.

For the calculation of the February temperature's time trend, we have also divided it into three periods. In the first period from 1915 to 1964, there is absolutely no time trend that reaches a significant level. In the other two periods, although there are regions with warming that passes the significance level, the warming areas between the two periods do not overlap, and their warming magnitude and range are both smaller than in January.

Wednesday, August 9, 2023

Manitoba: Analysis for the Month of January

The weather observation data from the Manitoba meteorological station sites have been used as the input, and the local spatial distribution of temperature data for over a hundred continuous months has been obtained through ANUSPLIN interpolation. This is the foundation of our data analysis in this series.

This is the first part of the series, a January analysis, with a time range from 1915 to 2023. Temperatures are represented with deep red for higher temperatures and deep blue for lower temperatures. As shown in the following chart, over each decade, the frequency of warmer January temperatures has notably increased compared to cooler ones.

The temperature in January has been continuous for a total of 109 years, and our research divides them into three periods to calculate and identify their time series trends. From the calculated results, there is absolutely no significant time series trend in the study area between 1915 and 1964. In the other two periods, as shown in the following chart, the vast majority of areas show a significant warming trend. It was also found that the warming trend between 1965 and 2023 is noticeably higher than between 1915 and 2023.

Sunday, August 6, 2023

ANUSPLIN: Configuration of Interpolation Range

The original elevation has a spatial resolution of 0.0041666665°, with a total of 3124 columns and 2641 rows, as shown in the figure below:

Observing the above figure, we can see that the longitude range is from -102.0083363186 to -88.9916701726, and the latitude range is from 48.99791691515 to 60.00208314165. If the cmd file of the lapgrd process directly uses the mentioned range as the interpolation range, it will result in an error. The reason for this error is that:

1 2	-102.0083363186+31240.0041666665=-88.9917 48.99791691515+26410.0041666665=60.0021

So, the range in the cmd file of lapgrd should be represented as follows:

1 2	-102.0083363186 -88.9917 0.0041666665 48.99791691515 60.0021 0.0041666665

If the error message 'GETLIM - LIMITS EXCEED SURFACE COORDINATE LIMITS' is encountered, the root cause of the issue may be that the range set in the 'splina' command is smaller than the actual DEM (Digital Elevation Model) range. To resolve this, simply adjust the range in the 'splina' command to be slightly larger to cover the entire extent of the DEM. This should allow the interpolation process to proceed without encountering the error.

Share another successful case study:

Lapgrd:

1 2	0.29531831E+06 0.17793583E+07 30 0.24045270E+07 0.40906770E+07 30

When the spatial extent of the lapgrd's cmd file is difficult to determine, consider specifying the spatial range precisely in the splina's cmd file, and use the .sur file obtained from splina's computation to determine the spatial extent of lapgrd.

当插值温度（或者其他情况需要表面协变量的）时，需要同时输出多个表面结果，需要注意此时的DEM文件在第二步的参数中需要多次设置，设置数量要与输出数量一致。

TEM11.sur
0
1
TEM11.cov
0
1
1
8871611.457651 10396111.457651 250
2
5131802.4564656 6137302.4564656 250
0
2
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
dem1.txt
2
-99.0
tem200001.grd
tem200002.grd
tem200003.grd
tem200004.grd
tem200005.grd
tem200006.grd
tem200007.grd
tem200008.grd
tem200009.grd
tem200010.grd
tem200011.grd
tem200012.grd
(100f10.3)

不建议批量插值，可能会出现意想不到的错误，推荐是一次插值一张而已。

Wednesday, July 27, 2016

Matlab+ANUSPLIN: ANUSPLIN Apps

List of Applications

[1] 太阳辐射下载代码，自URL批量下载；

[2] 太阳辐射数据转至XLS；

[3] 温度、降水数据由原始状态输出为XLS文件，注意这里未做有效性检验，早年数据质量不高，需特别注意.

[4] 温度ANU插值代码.

[5] 降水ANU插值代码.

[6] 太阳辐射ANU插值代码（无降水协变量）.

[7] 月值气象数据累积（或平均）转至年值数据，须包含mask文件.

[8] 多年（多期）气象数据求取平均值.

[9] 气象数据计算距平，须预先有均值数据.

[10] 太阳辐射在两种单位之间自由切换.

[11] 依据站点编码挑选坐标.

[12] 气象站属性（2019.04.13，TotalHub.xlsx）.

[13] 温度和降水插值（输出标准误差）.

Wednesday, June 15, 2016

Matlab+GDAL: ANUSPLIN数据处理套装

Introduction

ANUSPLIN插值过程需要做一些数据格式方面的转换，这里就是为它专门制作一个M代码组合套装，调用GDAL函数。

MAIN_TAG_NODATA_GTIFF.m给原Geotiff文件加上指示Nodata的Tag，方便后续转出为Arc/Grid时有NODATA_value这一配置。

MAIN_SETTING_PRJ_FILE.m带给ANUSPLIN输出的单独.grd文件加上投影信息（当然，他们的投影是一致的），方便随后转出为GeoTiff文件。

MAIN_CONVERT_FORMAT_SOLUTION.m是文件格式转换的解决方案，注意以下转换参数所对应的情况，及时调整，GDAL文件格式列表见参考文献[1]。

1
2
3

tarfot='GTiff';
ext='.tif';
files=dir(fullfile(SouDir, '*.grd'));

File.rar.注意查看.prj文件，这里默认参数是WGS_1984_UTM_Zone_50N.

References

[1] GDAL Raster Formats.

Tuesday, March 15, 2016

Design: ANUSPLIN空间插值植被指数

Fig. 1

File.rar. Forged NDV Code.

Fig. 2

Monday, December 14, 2015

Matlab: Conversion from GeoTiff Format to Arc/Info Grid

Summary

我在之前的帖子已经讨论过Arc/Info Grid或IDRISI Grid转换为GeoTiff格式文件的方法。此处，本帖讨论Geotiff文件转换为Arc/Info Grid文件的代码。

这里的方法只考虑到了像元经向、纬向空间分辨率一致的情况，对于这二者不一致的情况还没有做讨论，留待以后再考虑。

这一方法可以适用于任意投影的Geotiff文件。

Gtiff2ArcGrid.rar.

第二版.

Friday, November 20, 2015

Matlab: Converison from Arc/Info Grid to GeoTiff Format

Summary

笔记：记录ANUSPLIN LAPGRD输出.grd文件转换为通用文件Geotiff之方法。

Note

The codes in this case is to be oriented to the files associated with Lat/Long coordinate system only. For others, I would code in the future probably.

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

clear;
clc;

% Source Directory
SouDir='./Input';
% Destination Directory
DesDir='./Output';

% files
files=dir([SouDir, '/*.grd']);

for ii=1:length(files)
    filename=files(ii).name;
    filepath=[SouDir, '/', filename];
    
    % Output file fullpath
    outname=[filename(1:end-3), 'tif'];
    outpath=[DesDir, '/', outname];
    
    % Conversion
    Gonvert2GTiff(filepath, outpath);
    
    disp([num2str(ii), ':', outpath]);
end


disp('***************************************');

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

function Gonvert2GTiff(grdpath, gtiffpath);

    % Get the properties of this .grd file
    GeoInfo=GetProperty(grdpath);
    
    % Get the values of this .grd file;
    MatVal=GetValues(grdpath);
    MatVal=reshape(MatVal, GeoInfo.NCOLS, GeoInfo.NROWS);
    MatVal=MatVal';
    
    
    % Spatial extent
    stLong=GeoInfo.XLLCORNER;
    edLong=stLong+GeoInfo.CELLSIZE*GeoInfo.NCOLS;
    stLat=GeoInfo.YLLCORNER;
    edLat=stLat+GeoInfo.CELLSIZE*GeoInfo.NROWS;
    
    % Write a Data Referenced to Geographic Coordinates
    R=georasterref('RasterSize', [GeoInfo.NROWS, GeoInfo.NCOLS], 'Latlim', [stLat, edLat], 'Lonlim', [stLong, edLong], 'ColumnsStartFrom', 'north');
    
    geotiffwrite(gtiffpath, MatVal, R);



end

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Get values of a .grd file

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

function Output=GetValues(grdpath)
    
    % Read values
    fid=fopen(grdpath);
    
    tline=fgets(fid);
    count=1;
    Output=[];
    while ischar(tline)
        
        if count>6
            Vals=str2num(tline);
            Output=[Output; Vals'];
        end       
        tline=fgets(fid);
        count=count+1;
    end
    
     fclose(fid);
    fclose('all');


end

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy


function GeoInfo=GetProperty(grdpath)
    

    fid=fopen(grdpath);
    tline=fgets(fid);
    for ii=1:6
        [property, value]=strread(tline, '%s %f');
        str_exe=['GeoInfo.', upper(property{1}), '=', num2str(value), ';'];
        eval(str_exe);
        tline=fgets(fid);
    end
    
    fclose(fid);
    fclose('all');
    
end

References

[1] 实现气象数据的空间分布.

Monday, November 16, 2015

Matlab: Conversion from Arc/Info Grid to GeoTiff Format Program

Summary

记录ANUSPLIN插值输出结果在投影坐标系下的ASCII文件转换为GeoTiff格式文件。在前面的帖文中，我已经介绍Arc/Info Grid（Geographic Coordinate System）转换为Geotiff代码，此处追加Arc/Info Grid（Projected Coordinate System）转换为Geotiff代码。

Note

代码要求必须有一个Geotiff文件，它包含待转换Arc/Info Grid文件的全部地理信息，这些信息将在代码中赋给转换的矩阵进而输出为Geotiff文件。

Fig. 1

% Created by LI Xu
% Version 1.0
% October 26, 2015

% Description
% Main program for conversion to geotiff format 
% from Arc/Info Grid

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

clear;
clc;

% Mask
maskpath='mask.tif';
[~, geo]=geotiffread(maskpath);
info=geotiffinfo(maskpath);

% Source Directory
SouDir='./input';
% Destination Directory
DesDir='./output';

% All files
files=dir([SouDir, '/*.grd']);

for ii=1:length(files)
    filename=files(ii).name;
    % If exist
    otname=strsplit(filename, '.');
    otname=otname{1};
    otpath=[DesDir, '/', otname, '.tif'];
    
    if exist(otpath, 'file')
        disp([num2str(ii), ': ', otname, '.tif']);
        continue;
    end
    
    
    filepath=[SouDir, '/', filename];
    values=GetValues(filepath);
    geoattrs=GetProperty(filepath);
    inmat=reshape(values, geoattrs.NCOLS, geoattrs.NROWS);
    inmat=inmat';
    
    % Export
    
    geotiffwrite(otpath, inmat, geo, 'GeoKeyDirectoryTag', info.GeoTIFFTags.GeoKeyDirectoryTag);
    disp([num2str(ii), ': ', otname, '.tif']);
end

disp('****************************************************');

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Convert a raster into GeoTiff format from .grd files,
% generated by LAPGRD in ANUSPLIN.


% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy


function GeoInfo=GetProperty(grdpath)
    

    fid=fopen(grdpath);
    tline=fgets(fid);
    for ii=1:6
        [property, value]=strread(tline, '%s %f');
        str_exe=['GeoInfo.', upper(property{1}), '=', num2str(value), ';'];
        eval(str_exe);
        tline=fgets(fid);
    end
    
    fclose(fid);
    fclose('all');
    
end

% Created by LI Xu
% Version 1.0
% 10 September, 2015

% Description:
% Get values of a .grd file

% If you have any question about this code,
% please do not hesitate to contact me via E-mail: 
% jeremy456@163.com

% Blog:
% http://blog.sciencenet.cn/u/lixujeremy

function Output=GetValues(grdpath)
    
    % Read values
    fid=fopen(grdpath);
    
    tline=fgets(fid);
    count=1;
    Output=[];
    while ischar(tline)
        
        if count>6
            Vals=str2num(tline);
            Output=[Output; Vals'];
        end       
        tline=fgets(fid);
        count=count+1;
    end
    
    fclose(fid);
    fclose('all');

end

Attachment.

LI Xu's World

Sunday, October 15, 2023

Batch interpolation of meteorological variables with ANUSPLIN

References

Saturday, August 12, 2023

Manitoba: Analysis for the Month of February

Wednesday, August 9, 2023

Manitoba: Analysis for the Month of January

Sunday, August 6, 2023

ANUSPLIN: Configuration of Interpolation Range

Wednesday, July 27, 2016

Matlab+ANUSPLIN: ANUSPLIN Apps

List of Applications

Wednesday, June 15, 2016

Matlab+GDAL: ANUSPLIN数据处理套装

Introduction

References

Tuesday, March 15, 2016

Design: ANUSPLIN空间插值植被指数

Monday, December 14, 2015

Matlab: Conversion from GeoTiff Format to Arc/Info Grid

Summary

Friday, November 20, 2015

Matlab: Converison from Arc/Info Grid to GeoTiff Format

Summary

Note

References

Monday, November 16, 2015

Matlab: Conversion from Arc/Info Grid to GeoTiff Format Program

Summary

Note

Labels

Blog Archive