[ncl-talk] 回复： How to get the regional data of SMAP L4 *(HDF5 format)
穆二
dg1225033 at smail.nju.edu.cn
Tue Dec 19 19:15:00 MST 2017
Dear Dennis Shea,
It is really very nice of you to pay attention to the uploaded ncl script for extracting subset data of SMAP L4, your professional and detailed analysis and reply helped me a lot again, your analysis deepened my understanding for the missing value in NCL, and this is a basic and very important concept when writing script, especially for calculating.
Thanks a lot for your patient guidance and help. Greatly appreciate for your attention and time!
Christmas is fast approaching, I wish you and other members in ncl-talk Merry Christmas!
Cheers!^-^
Sincerely,
Yi Lu
------------------ 原始邮件 ------------------
发件人: "Dennis Shea";<shea at ucar.edu>;
发送时间: 2017年12月18日(星期一) 晚上11:09
收件人: "易路"<dg1225033 at smail.nju.edu.cn>;
抄送: "Ncl-talk"<ncl-talk at ucar.edu>;
主题: Re: [ncl-talk] How to get the regional data of SMAP L4 *(HDF5 format)
You are welcome.
It is best not to send direct emails.
Please send all emails to ncl-talk at ucar.edu (cc'd here)
==============
The original script had the comment:
;---Manually add a _FillValue to lat/lon; not sure why
;---_FillValue not associated with the variable
lat2d at _FillValue = -9999.0
lon2d at _FillValue = -9999.0
;printMinMax(lat2d, 0)
;printMinMax(lon2d, 0)
The comment should have been:
;---Manually add an _FillValue attribute to lat/lon
;---Not sure why but there is no _FillValue associated with the
;---'latitude' and 'longitude' variables on the file
Variable: lat2d <===== printVarSummary(lat2d)
Type: float
Total Size: 1565536 bytes
391384 values
Number of Dimensions: 2
Dimensions and sizes: [DIM_000 | 406] x [DIM_001 | 964]
Coordinates:
Number Of Attributes: 2
long_name : Latitude of the center of the Earth based grid cell.
units : degrees_north
===
Note: There is no _FillValue attribute associated with lat2d ('latitude' on the file).
However, if you use 'printMinMax(lat2d,0)', you will see a minimum value of -9999.
This is why ncl-talk frequently states:
The most important rule in data processing is: Look at your data
Using
printMinMax(lat2d,0) yields
(0) Latitude of the center of the Earth based grid cell. (degrees_north) : min=-9999 max=83.632
Obviously, -9999.0 is not a valid latitude.
'We' had to manually attach the _FillValue in the script to assure correct results.
Cheers
D
On Sun, Dec 17, 2017 at 2:56 AM, 易路 <dg1225033 at smail.nju.edu.cn> wrote:
Dear Dennis Shea,
Thanks very much for your such professional test to check whether the grids of SMAP L4 is rectilinear. This test result that the grids of SMAP L4 belongs to rectilinear solved my confusion why I could successfully extracte the regional data by coordinate subscripting once the soil moisture variable was given the coordinate.
Looking into the matrix of la2d and lon2d of the SMAP L4, I found that although they are two dimensional, but acturally they can be recorded only in one dimension, since each row of lon2d is the same, and each column of lat2d is the same. So, I wrote the script as follows to cordinate extract the regional data, and it worked.
lat2d = f->cell_lat
lon2d = f->cell_lon
soilMoisture = f->sm_rootzone
soilMoisture!1="lat"
soilMoisture!2="lon"
soilMoisture&lon=lon2d(0,:)
soilMoisture&lat=lat2d(:,0)
soilMoisture_region = soilMoisture({maxLatBoundary:minLatBoundary},{minLonBoundary:maxLonBoundary})
It's very nice of you to correct my understanding about the relation ship between lat/lon dimension and the grid structure. Two dimensions are 'almost always' but not absolutely curvilinear. This is very important for me!
Greatly appreciate for your precious attention, time and professional help.
with best regards!
Sincerely,
Yi Lu
易路
南大邮件系统/学生/博士生/12级博士生
南京市汉口路22号
------------------ Original ------------------
From: "Dennis Shea"<shea at ucar.edu>;
Date: Sat, Dec 16, 2017 01:03 AM
To: "Mary Haley"<haley at ucar.edu>;
Cc: "易路"<dg1225033 at smail.nju.edu.cn>; "ncl-talk"<ncl-talk at ucar.edu>;
Subject: Re: [ncl-talk] How to get the regional data of SMAP L4 *(HDF5 format)
What Mary wrote is generally correct ... Specifically:
========================
You are trying to extract the region from your HDF5 data by using coordinate subscripting.
Coordinate subscripting can only be done on rectilinear data, which is data represented by one-dimensional lat/lon arrays that are the same length as their corresponding dimension.
Your data is likely curvilinear, which means your data is represented by two-dimensional lat/lon arrays. I say this because in your code it looks like you are extracting 2D lat/lon arrays:
lat2d = f->cell_lat
lon2d = f->cell_lon
nSmapLat = dimsizes(lat2d(:,0))
nSmapLon = dimsizes(lon2d(0,:))
===========================
Generally, when people see two dimensional latitude & longitude arrays ( cell_lat[*][*] and cell_lon[*][*] ), the grid structure is 'almost always' curvilinear. However, I noted that the file is a L4 file (SMAP_L4). NASA Level-4 files are 'almost always' on a 'nice' grid.
I used an undocumented and unsupported function of mine to test if the grid is a rectilinear grid despite the two-dimensional latitude & longitude arrays. The function returns True if it is rectilinear and False otherwise.
%> ncl tst.smap_rectilinear.ncl
The grid is rectilinear. Hence, you can use the approach you originally used. (It did need a few improvements :-) )
%> ncl smap_regional.ncl_YiLu
Yes ... there are a fe locations over the water but if you use a higher resolution map that will not be the case.
Good luck
On Thu, Dec 14, 2017 at 9:21 AM, Mary Haley <haley at ucar.edu> wrote:
Dear Yi Lu,
You are trying to extract the region from your HDF5 data by using coordinate subscripting.
Coordinate subscripting can only be done on rectilinear data, which is data represented by one-dimensional lat/lon arrays that are the same length as their corresponding dimension.
Your data is likely curvilinear, which means your data is represented by two-dimensional lat/lon arrays. I say this because in your code it looks like you are extracting 2D lat/lon arrays:
lat2d = f->cell_lat
lon2d = f->cell_lon
nSmapLat = dimsizes(lat2d(:,0))
nSmapLon = dimsizes(lon2d(0,:))
Given this information, you cannot use coordinate subscripting to subscript your data. Instead, you either need to use the function getind_latlon2d or region_ind to get the indexes of the area of interest. For example:
lat_reg = (/minLatBoundary,maxLatBoundary/)
lon_reg = (/minLonBoundary,maxLonBoundary/)
ij = getind_latlon2d(lat2d,lon2d,lat_reg,lon_reg)
sm_regional = sm(ij(0,0):ij(1,0),ij(0,1):ij(1,1))
Before plotting, you need to subset your lat2d / lon2d arrays using the same indexing:
res at sfXArray = lon2d(ij(0,0):ij(1,0),ij(0,1):ij(1,1))
res at sfYArray = lat2d(ij(0,0):ij(1,0),ij(0,1):ij(1,1))
Please see our "Subsetting / extracting data based on lat / lon values" example page:
http://www.ncl.ucar.edu/Applications/latlon_subset.shtml
and in particular, look at examples latlon_subset_2.ncl and latlon_subset_3.ncl, which show the difference between getind_latlon2d and region_ind.
Also, to better understand the different types of data (rectilinear, curvilinear, unstructured) and how to plot them, see:
http://www.ncl.ucar.edu/Applications/plot_data_on_map.shtml
Good luck,
--Mary
On Wed, Dec 13, 2017 at 7:05 PM, 易路 <dg1225033 at smail.nju.edu.cn> wrote:
Hi, all
Sorry to disturb you for technical help.
I am trying to extract the regional data from the soil product SMAP L4 which is stored in the format of HDF5. I wrote the ncl script taking referecnce to the example listed on the web of http://www.ncl.ucar.edu/Applications/HDF.shtml. But the extracted data were not in the demond region. So would you please help me to find the faults in my ncl script? The main part of this script is as follow, and the detail information can be downloaded from the accessory.
********main part of the script *********
nFiles = dimsizes(iFiles)
f = addfile(iFiles(0), "r")
lat2d = f->cell_lat
lon2d = f->cell_lon
nSmapLat = dimsizes(lat2d(:,0))
nSmapLon = dimsizes(lon2d(0,:))
minSmapLat = min(lat2d)
maxSmapLat = max(lat2d)
minSmapLon = min(lon2d)
maxSmapLon = max(lon2d)
lat=fspan(minSmapLat, maxSmapLat, nSmapLat)
lat!0="lat"
lat&lat=lat
lat at units="degrees_north"
lat at long_name="Latitude"
lon=fspan(minSmapLon,maxSmapLon, nSmapLon)
lon!0="lon"
lon&lon=lon
lon at units="degrees_east"
lon at long_name="Longitude"
lat2d!0="lat"
lat2d!1="lon"
lat2d&lon=lon
lat2d&lat=lat
lon2d!0="lat"
lon2d!1="lon"
lon2d&lon=lon
lon2d&lat=lat
minLatBoundary = 25
maxLatBoundary = 35
minLonBoundary = 110
maxLonBoundary = 120
sm = f->sm_rootzone
printVarSummary(sm)
sm!0="lat"
sm!1="lon"
sm&lon=lon
sm&lat=lat
printVarSummary(sm)
sm_regional = sm({minLatBoundary:maxLatBoundary},{minLonBoundary:maxLonBoundary})
***********************************
Thanks a lot for your precious attentions and time!
Yi Lu
易路
南大邮件系统/学生/博士生/12级博士生
南京市汉口路22号
从腾讯企业邮箱发来的超大附件
SMAP_L4_SM_gph_20150815T013000_Vv3030_001.h5 (136.33M, 2018年01月12日 11:27 到期)进入下载页面：http://exmail.qq.com/cgi-bin/ftnExs_download?k=523465623d8c4a92182f6ffe1038075719575c5a030c065706190304015b1d51000c504f0e0e55561b075257010c0452050d070051381e65657924326874043a65793a0547506f57060550520f0905310605565207086f334007555107670055181a4b4c5f0d3058&t=exs_ftn_download&code=64eb780e&s=email
_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
_______________________________________________
ncl-talk mailing list
ncl-talk at ucar.edu
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20171220/e065c946/attachment.html>
More information about the ncl-talk
mailing list