Transforming the default projection system

Hi,

Good afternoon. I have a NetCDF file of daily PM2.5 concentration. The file is projected in Lambert Conformal Conic projection, and now I want to project it into WGS84. I am wondering is there any post-processing tool/ioapi script exists that allows converting a NetCDF file into another projection. Or, some sort of script, which can convert the latitude and longitude of a given projection (csv, netcdf or other format) to others.

I would highly appreciate it if you could share your knowledge to help me came out of this issue.

Thanks,
Sadia

See:

You don’t have things quite right in your question: WGS84 is a geodetic spheroid, a particular choice of representing the Earth’s surface as an oblate ellipsoid; WGS-84 os one of these and the WRF-Sphere is another.
Lambert Conformal Conic is a map projection, a scheme for putting 2-D Cartesian X-Y coordinates on a (portion of) the Earth’s surface (the transformation equations assume that you have already chosen a geodetic spheroid, btw. Different spheroids may well give differences in coordinates as large as 25KM at CONUS scales.)

Grids are defined in terms of a map projection: given the projection, specify the lower-left corner coordinates, the cell-sizes DX in the X direction and DY in the Y direction, and the numbers NCOLS of grid-columns and NROWS of grid-rows in order to specify a grid.

Various of the m3tools programs can be used to transform gridded data from one grid to another (each of which has potentially its own map-projection. Program m3cple does linear interpolation, which is appropriate for going from a coarse (large-cells) grid to a find (small-cells) grid. The pair of programs mtxcalc+mtxcple does re-sampling and re-aggregation, appropriate for going from a finer grid to a coarser grid (or grids of about the same cell-size); it also works for going from coarser to finer.

BEWARE that coordinate-transformations for emissions data have a piece of non-obvious trickiness: as per the usual conventions, emissions data is not given in terms of standard MKS units, but instead has an extra (implicit) per-grid-cell qualifier to the units: units recorded as Kg/sec really mean Kg/sec/grid-cell, so in doing transforms of emissions data you need also to adjust for the ratio of the output-grid grid-cell area to the input-grid grid-cell area. Program mtxcalc has an option for doing this adjustment, so you need to use mtxcalc+mtxcple for re-gridding emissions data.

1 Like

Hi,

Thank you so much for your reply. It was really very helpful post. Thanks for your time for pointing to some helpful resources and making the concept clear to me.

Regards,
Sadia

Hi Carlie,

Sorry for bothering you again. I was reading the link(https://cjcoats.github.io/ioapi/GRIDS.html) you shared. I have a question. If now I understand correctly, " Grids are defined in terms of a map projection", for my case it is “Lambert Conformal Conic”. But I didn’t understand the significance of the parameter “GCTP projection 4”. Also, how can I know which “geodetic spheroid” has been chosen? Is there any way to know that? I would appreciate your help in this regard. Thanks again for clarifying my concept on projection.

Regards,
Sadia

OK. let’s take it from the top.

The I/O API coordinate system treatment is based upon the USGS Geographic Coordinate Transform Package (GCTP), which at the time we started the project was the easiest high-quality publicly-available geographic package for us to use. Its original manual is at https://cjcoats.github.io/ioapi/GCTP.pdf and many of the I/O API coordinate concepts and parameters come from it.

The Earth is approximately an ellipsoid, in particular, because of the force of the Earth’s rotation it is a “flattened” spheroid – the rotation causes the diameter to “stretch” in the plane of the Equator (so that the diameter across the Equator is greater than the diameter pole-to-pole. A spheroid is such a description of the Earth; there have been many such descriptions proposed and used; the most common, together with their GCTP code-numbers are

  • Clarke 1866: GCTP uses 0 as its index-number
  • GRS 1980: 8
  • WGS 84: 12
  • perfect sphere: 19 (in which case we have to specify the radius elsewhere)
    Generally the more recent GRS 80 or WGS 84 ellipsoids are what GIS people use; meteorologists tend to use perfect spheres since doing proper equations of motion on a more-general ellipsoid is quite “hairy” mathematics. I have seen some older data given in terms of Clarke 1866…

A map projection is a way of putting Cartesian X-Y coordinates on ellipsoids (and so needs a ellipsoid-specification to complete its definition). A Lambert Conformal Conic, for example (I/O API projection code 2, GCTP projection code 4), has the geometric algorithm:

  • Form a cone that intersects the ellipsoid at two latitudes alpha and beta (and has central axis the same as the rotational N-S axis of the Earth).
  • Project a geometric ray from the center of the Earth through both its surface and the surface of the cone. Map the Earth-surface intersection-point of that ray to the cone’s intersection point. If you do the arithmetic for this you’ll get into some rather messy trigonometry :slight_smile:
  • Slice the cone diametrically opposite “central longitude” gamma and unroll it flat, with the longitude-gamma’s image on the cone in the N-S, or Y direction.
  • Find the point on the unrolled-cone where “origin-point” latitude-longitude (LATORIG,LONORIG) winds up, and start measuring X and Y coordinates from there – Y parallel to the longitude-gamma’s line. Usually one chooses LONORIG=gamma, and LATORIG somewhere between alpha and beta (but not always;-( )
    Because of how this works, choosing different ellipsoids to represent the Earth will give somewhat-different X-Y coordinate systems; if you use Lambert projections with alpha=30, *beta=*60, and gamma=90, LONORIG=-90 (note the minus-sign, for longitudes west of the Prime Meridian LON=0), and LATORIG=45, then the Cartesian origin X=Y=0 somewhere in Illinois no matter which ellipsoid you use, but the coordinates for Portland Maine will be different by about 20 KM between choosing a spherical ellipsoid and the WGS84 ellipsoid.

Note that because of round-off problems in all the map-projection arithmetic, it is best to do everything in REAL8 (DOUBLE PRECISION) rather than in REAL4, which doesn’t have enough significant digits to be sure of giving the “right” answer…

And the I/O API’s coordinate transform routines generally do this work for you, given descriptions you don’t need to worry much about the messy mathematics.

The worst problem you’re likely to encounter is transforming back and forth between grid-cell columns and rows and Cartesian X-Y coordinates:
Suppose you have a grid with lower-left corner (X0,Y0) and cell-size (DX,DY), and a point (X,Y) with X≥X0, Y≥Y0. Then the column-
and row-numbers (C,R) for the point are
C=1+INT((X-X0)/DX), R=1+INT((Y-Y0)/DY).
Beware that computer-languages “convert to integer” operators tend to “do the wrong thing” from a mathematical point of view: INT(-0.5) and INT(0.5) are both 0 ;-( so you have to do explicit checks for X≥X0, Y≥Y0 and say “the point is off the grid” if either one fails.

Thank you so much for your step by step explanation. I appreciate your time to make these things clear to me. Though I am still struggling with my issues, many things are clear to me now.

So, in summary, as my GCTP projection code is 4 according to https://cjcoats.github.io/ioapi/GRIDS.html, I can say that the sphere used for projection is “4. International 1909” [https://cjcoats.github.io/ioapi/GCTP.pdf].

Since I need to know the latitude and longitude associated with WGS 84 ellipsoid (12. WGS 84) from my GRIDDESC file. I was trying the following programs:

…GRIDDESC…
‘LCCPROJ’ 2 33.000 45.000 -97.000 -97.000 40.000
‘4NC3’ ‘LCCPROJ’ 960000.000 -660000.000 4000.000 4000.000 237 153 1
…GRIDDESC…

-->setenv GRIDDESC /home/safrin/GRIDDESC
-->setenv IOAPI_ISPH 12
-->setenv GRIDDED /home/safrin/gridded_file
-->setenv BOUNDARY "NONE"
-->latlon GRIDDESC GRIDDED BOUNDARY >& summary

But, I am getting the following error.
“latlon: error while loading shared libraries: libnetcdf.so.11: cannot open shared object file: No such file or directory”

If I am on the right track, do I need to install another version of NetCDF?

I am feeling really very sorry for bothering you a lot.

Thanks,
Sadia

Where did your copy of the program latlon come from?

Did you build it yourself (or did someone build it for you) according to the directions in https://cjcoats.github.io/ioapi/AVAIL.html#build?

Those instructions told you how to build a statically linked executable, that would not need that shared-object file.

No, I have not built the latlon (https://cjcoats.github.io/ioapi/LATLON.html). Thanks for the link. I will follow the instructions.

Regards,

Sadia