Re: re-converting 1D to 3D array.

Sho,

My mistake, I did not notice the 2-D coordinates before. There is
an efficient method to handle subsetting for 2-D coordinates. Use
the region_ind function to get the bounding x and y subscripts for
the region of interest:

http://www.ncl.ucar.edu/Document/Functions/Contributed/region_ind.shtml

Here is an untested upgrade to my first suggestion. This gets the
bounding subscript values, then reads a 3-D subset data array of the
desired dimensions.

   latS = 37
   latN = 47
   lonW = 269
   lonE = 279

   lat = fin->lat ; read 2-D coordinates
   lon = fin->lon

   ji = region_ind (lat, lon, latS, latN, lonW, lonE)
   j1 = ji(0) ; lat start ; get subset indices
   j2 = ji(1) ; lat last
   i1 = ji(2) ; lon start
   i2 = ji(3) ; lon last

   prc = fin->pr(:,{i1:i2},{j1:j2}) ; read 3-D subset

   lats_region = lat({i1:i2},{j1:j2}) ; make matching subset
   lons_region = lon({i1:i2},{j1:j2}) ; of 2-D coordinates

   printVarSummary (prc) ; check results before proceeding
   print (lats_region)
   print (lons_region)
   exit

There may be be a few extra grid points along the edges where the
subset data grid is out of alignment with a perfect rectangle as
defined with lat/lon boundaries. You can then set these edge points
to missing values if needed, but it is simpler for many applications
to just leave them alone. I hope this is the right solution for you.

--Dave A.
NOAA/PSD/CIRES

Sho Kawazoe wrote:
> Hi,
>
> And thanks for the quick reply. However, based on your answer,
> this assumes that the pr variable in the NetCDF file is originally in (time,
> lat, lon) dimensions correct? It is in (time, xc, yc) coordinate, so is
> there an maneuver to get around this?
>
> Sho
>
> On Tue, Mar 9, 2010 at 6:46 PM, Dave Allured <dave.allured@noaa.gov> wrote:
>
>> Sho Kawazoe,
>>
>> This is the second statement in your program:
>>
>> prc = fin->pr
>>
>> Move up the four lines for lat1, lat2, lon1, lon2; and try this instead:
>>
>> prc = fin->pr(:,{lon1:lon2},{lat1:lat2})
>>
>> printVarSummary(prc)
>> print(prc&lon)
>> print(prc&lat)
>> exit
>>
>> If latitudes in the file go from high to low, then reverse the order
>> of lat1 and lat2 above. Same for longitudes.
>>
>> Somehow you started making this reader more complicated than it
>> should be. This should work better. We do not need to discuss
>> conversion of dimensions unless your problem is more complicated
>> than it looks.
>>
>> --Dave A.
>> NOAA/PSD/CIRES
>>
>> Sho Kawazoe wrote:
>>> NCL users,
>>>
>>> I'm still very raw in my knowledge with ncl, and was wondering if anyone
>>> can help me out.
>>>
>>> I'm currently using precipitation data with MM5 (from EGS), and trying to
>>> mask a specific lat and lon point, with the ultimate purpose of
>> outputting
>>> data within that specific grid, and truncating the data to where when I
>> end
>>> up concatenating other files, I am well below the 2gb threshold.
>>>
>>> Steps leading to retrieving the desiredvalues in 1D looks fine, and the
>>> output file is well below the original input file size. However, when I
>>> reconvert the 1D to 3D, the file is back to the same size as the input
>> file,
>>> and it looks like missing values are inputted into the output file, which
>>> doesn't make the file any smaller.
>>>
>>> Is there a way to reconvert the 1D array to 3D array and keep only the
>>> precip values we want, without the missing values(keep in mind that the
>>> associated metadata needs to be there as well)? A simple onedtond doesn't
>>> look like the solution, and I'm misunderstanding the coding.
>>>
>>> The script is as follows. Any help is appreciated.
>>>
>>> ;*****************************************************
>>> ; Read original NetCDF file.
>>> ;*****************************************************
>>>
>>> fin = addfile ("pr_MM5I_1979010103.nc","r")
>>>
>>> ;******************************************************
>>> ; Upload netCDF file variables, set dimensions
>>> ;******************************************************
>>>
>>> prc = fin->pr ; obtain pr(time,xc,yc) values. Dimension is
>> in
>>> 3D
>>> lat = fin->lat
>>> lon = fin->lon
>>> time = fin->time
>>> Lambert_Conformal = fin->Lambert_Conformal
>>>
>>> dim = dimsizes(prc)
>>> ntime = dim(0)
>>> nlat = dim(1)
>>> nlon = dim(2)
>>>
>>> latlonmask2 = prc(:,:,:) ; assign dim size to masking var
>>> latlonmask2 = 0 ; initialization
>>>
>>> ;**************************************************************
>>> ; Indicate coordiate you want to analyze. Make sure you
>>> ; check the NetCDF file for correct lat lon specification!!
>>> ;
>>> ; Mask through time, xc, and yc for entire precip variable to
>>> ; analyze specifc regions
>>> ;**************************************************************
>>>
>>> lat1 = 37
>>> lat2 = 47
>>> lon1 = 269
>>> lon2 = 279
>>>
>>> do k = 1, ntime-1
>>> do j = 1, nlat-1
>>> do i = 1, nlon-1
>>> if (lat(j,i).ge.lat1 .and. lat(j,i).lt.lat2 .and. lon(j,i).gt.lon1
>>> .and. lon(j,i).le.lon2) then
>>> latlonmask2(k,j,i) = 1
>>> end if
>>> end do
>>> end do
>>> end do
>>>
>>> desiredvalues = mask(prc, latlonmask2, 1)
>>>
>>> ;********************************************************************
>>> ; Convert the masked region to a 1d array, which allows us to then
>>> ; use the ncl function "ind" to keep only the data from the region
>>> ; we are interested in looking at. This also keeps all the
>>> ; metadata associated with each precip value.
>>> ;*******************************************************************
>>>
>>> desired1d = ndtooned (desiredvalues) ; convert to 1d array
>>> indices = ind (.not.ismissing(desired1d)) ; keep precip values that
>>> have no missing data.
>>> output1d = desired1d(indices) ; all non-missing values
>>> desiredvalues = onedtond(output1d,dimsizes(desiredvalues))
>>> delete (desired1d) ; clean up if you wish.
>>> delete (indices)
>>>
>>> ;******************************************************************
>>> ; Create a new nc output file
>>> ;******************************************************************
>>>
>>> system("/bin/rm -f region_output.nc") ; remove any
>> pre-existing
>>> file
>>> fout = addfile("region_output.nc","c")
>>> filedimdef(fout,"time",-1,True) ; Makes time an UNLIMITED
>>> dimension (as suggested by NCL)
>>>
>>> fout->prc = output1d
>>>
>>> fout->lat = lat
>>> fout->lon = lon
>>>
>>> end
>>>
>>>
>>> ------------------------------------------------------------------------
>>>
>>> _______________________________________________
>>> ncl-talk mailing list
>>> List instructions, subscriber options, unsubscribe:
>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>> _______________________________________________
>> ncl-talk mailing list
>> List instructions, subscriber options, unsubscribe:
>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>
>
>
>
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> ncl-talk mailing list
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
_______________________________________________
ncl-talk mailing list
List instructions, subscriber options, unsubscribe:
http://mailman.ucar.edu/mailman/listinfo/ncl-talk
Received on Wed Mar 10 13:50:10 2010