[ncl-talk] grib projection conversion

[Rick Brownrigg]

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
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20191204/c1a74de6/attachment.html>