[ncl-talk] grib projection conversion

[Micah Sklut]

I was able to find the wgrib2 tool, that has the new_grid option for
converting grib files projections. This works quite easy.
https://www.cpc.ncep.noaa.gov/products/wesley/wgrib2/new_grid.html

On Wed, Dec 4, 2019 at 1:35 PM Micah Sklut <micahs2005 at gmail.com> wrote:

> Thanks.  Still having trouble with this one. I get a return array in
> meters from the transform_coordinate function. It doesn't look like I can
> get output in degrees, so in order to convert my existing lat/lon grid into
> the new x/y coordinates in meters, I would have to convert my initial
> lat/lon coords to meters, and then I could use linint2?
>
>
> On Wed, Dec 4, 2019 at 12:25 PM Rick Brownrigg <brownrig at ucar.edu> wrote:
>
>> Yes, meters are the default units
>>
>> On Wednesday, December 4, 2019, Micah Sklut <micahs2005 at gmail.com> wrote:
>>
>>> Perhaps it's converting to meters rather than degrees. The poles (90,
>>> -90), come out as infinity, which makes sense for Mercator. I'll check out
>>> the proj documentation to see what output options there are.
>>>
>>> On Wed, Dec 4, 2019 at 11:56 AM Micah Sklut <micahs2005 at gmail.com>
>>> wrote:
>>>
>>>> Thanks. I'll look at ndtooned(), but I did get the function to work, it
>>>> just returns values that are way outside of the expected lat/lon range.
>>>>
>>>> Here's what my script looks like:
>>>>
>>>> data = "~/wc_data/data/"
>>>>
>>>> gfsFile = data + "gfs/gfs.t12z.pgrb2.0p25.f000.grb"
>>>>
>>>> f = addfile(gfsFile,"r");
>>>> lat = f->lat_0;
>>>> lon = f->lon_0;
>>>>
>>>> arraySize = dimsizes(lat) * dimsizes(lon)
>>>> x = new( (/ arraySize /), double)
>>>> y = new( (/ arraySize /), double)
>>>> z = new( (/ arraySize /), double)
>>>>
>>>> latSize = dimsizes(lat)
>>>> lonSize = dimsizes(lon)
>>>>
>>>> do j = 0, dimsizes(lat)-1
>>>> coordStart = j * lonSize
>>>> coordFinish = coordStart + lonSize - 1
>>>> x(coordStart:coordFinish) = lon(:)
>>>> y(coordStart:coordFinish) = lat(j)
>>>> z(coordStart:coordFinish) = 0
>>>> end do
>>>>
>>>> ;transform_coordinate(srcProj:string, dstProj:string, x[*]:double,
>>>> y[*]:double, z[*]:double)
>>>> transform_coordinate("+proj=lonlat +ellps=sphere", "+proj=merc", x, y,
>>>> z)
>>>>
>>>> print(x)
>>>>
>>>> On Wed, Dec 4, 2019 at 11:50 AM Rick Brownrigg <brownrig at ucar.edu>
>>>> wrote:
>>>>
>>>>> Yes, you'll need to flatten your arrays into 1D. See the ndtooned()
>>>>> and related functions for ways to do that.
>>>>>
>>>>> On Wednesday, December 4, 2019, Micah Sklut <micahs2005 at gmail.com>
>>>>> wrote:
>>>>>
>>>>>> Thanks Rick. I think this transform_coordinate function will take
>>>>>> care of my needs.
>>>>>>
>>>>>> In my case the gfs data has a latitude array of 1440, and longitude
>>>>>> array of 721.  So, does this mean, I need to pass transform_coordinate a 1d
>>>>>> array for X that shows latitude values for every coordinate? Meaning
>>>>>> 1440x721 values? And, subsequently for Y? And, I"ll pass a 1d array for Z
>>>>>> with all zeroes that have 1440x721 values.
>>>>>>
>>>>>> On Tue, Dec 3, 2019 at 5:13 PM Rick Brownrigg <brownrig at ucar.edu>
>>>>>> wrote:
>>>>>>
>>>>>>> Hi,
>>>>>>>
>>>>>>> NCL has an undocumented interface to the Proj4 cartographic
>>>>>>> projection library that could produce coordinates in a Mercator projection,
>>>>>>> given a list of lat/lons:
>>>>>>>
>>>>>>>    transform_coordinate(srcProj:string, dstProj:string, x[*]:double,
>>>>>>> y[*]:double, z[*]:double)
>>>>>>>
>>>>>>> Where srcProj/dstProj are the "proj4 strings" describing the source
>>>>>>> and destination projections. In your case, the srcProj would be something
>>>>>>> like:
>>>>>>>   "+proj=lonlat +ellps=sphere"     # I've encountered some instance
>>>>>>> where things fail without specifying the ellipsoid
>>>>>>>
>>>>>>> dest:
>>>>>>>
>>>>>>>   "+proj=merc"
>>>>>>>
>>>>>>> The x/y/z arrays are the same length; the z can be set to zeros if
>>>>>>> real data unavailable. The conversion happens "in place".  More on proj4
>>>>>>> here. Particularly if your area of interest is not global, there are likely
>>>>>>> other parameters to the Mercator you might want to add, like center of
>>>>>>> projection, etc:
>>>>>>>
>>>>>>>     https://proj.org/
>>>>>>>
>>>>>>> (Note, NCL's binding to proj4 is specifically the older v4.xxx
>>>>>>> libraries, not the newer v.6. series. The proj-strings should be compatible
>>>>>>> however)
>>>>>>>
>>>>>>> Finally, I don't understand the bit about starting with a 1440x721
>>>>>>> and ending up with a 1440x1440. It almost sounds like you'll need to regrid
>>>>>>> to that density in lat/lon space prior to projecting.
>>>>>>>
>>>>>>>
>>>>>>> Rick
>>>>>>>
>>>>>>>
>>>>>>> On Tue, Dec 3, 2019 at 3:01 PM Micah Sklut via ncl-talk <
>>>>>>> ncl-talk at ucar.edu> wrote:
>>>>>>>
>>>>>>>> Hi,
>>>>>>>>
>>>>>>>> I'm interested in the best way to convert lat/lon and corresponding
>>>>>>>> data values from GFS grib data into Mercator projection values. The use
>>>>>>>> case is not for mapping in NCL, but taking an input grib file and
>>>>>>>> converting values to an output grib file.
>>>>>>>>
>>>>>>>> The GFS data I have is 0.25x0.25, so 1440x721 values, so I believe
>>>>>>>> I'm looking for the Mercator projection to provide 1440x1440 values, where
>>>>>>>> the lat/lons are spaced according to the Mercator projection. Is there a
>>>>>>>> conversion function in NCL that will do this?
>>>>>>>>
>>>>>>>> Thanks!
>>>>>>>>
>>>>>>>> --
>>>>>>>> Micah Sklut
>>>>>>>>
>>>>>>>> _______________________________________________
>>>>>>>> ncl-talk mailing list
>>>>>>>> ncl-talk at ucar.edu
>>>>>>>> List instructions, subscriber options, unsubscribe:
>>>>>>>> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
>>>>>>>
>>>>>>>
>>>>>>
>>>>>> --
>>>>>> Micah Sklut
>>>>>>
>>>>>>
>>>>
>>>> --
>>>> Micah Sklut
>>>>
>>>>
>>>
>>> --
>>> Micah Sklut
>>>
>>>
>
> --
> Micah Sklut
>
>

-- 
Micah Sklut
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191205/4e994397/attachment.html>