Re: adding lat/lon variables from projected coordinate variables

Thank you very much, this worked perfectly.

I'll try to summarize my problem at the solutions as I expect others might
run into similar issues.

Problem: I converted projected raster data to NetCDF using GDAL. Only
projected x and y fields exist in the new NetCDF output by GDAL. Need
lat/lon coordinates (or corners at least) for display/analysis with NCL.

Solution:
(1) Use NCL to loop through the projected coordinates and export each pair
as a new line in an ascii file (see David's example code). Only need to
change slightly to export, for example: system("echo " +
sprinti("%7i",x(j)) + " " + sprinti("%7i",y(i)) + " >> exports/proj_yx.txt")

(2) Translate these projected coordinates to longitude/latitude. David's
example used inproj, which I believe may be deprecated. It appears that
'proj -i' does the same task. An alternative, which I found easier is to
use gdaltransform. Note that gdaltransform takes the X coordinate first,
then Y (need to be changed from Dave's original code if using
gdaltransform):

gdaltransform -s_srs EPSG:3412 -t_srs WGS84 < proj_xy.txt >
prox_xy_toLL.txt

This feeds the projected (in EPSG 3412 polar stereographic) coordinates in
XY format to gdaltransform and exports to another ascii in longitude
latitude. You'll need to either delete the third ascii column (Z) or
account for in in the next steps. Can delete it with 'awk': awk
'{$NF="";print}' proj_xy_toLL.txt > proj_xy_toLL_2.txt

(3) Read these new geographic coordinates with NCL, assign to 2d lon/lon
variables and write to NetCDF as David's sp.ncl script does. To use Dave's
code as is you'll also need to swap these geographic coordinates back to YX
format (lat-lon) using awk again: awk '{t=$1; $1=$2; $2=t;print;}'
proj_xy_toLL2.txt > proj_yx_toLL.txt

Thanks again!

Leif

On Tue, Feb 19, 2013 at 7:21 PM, David Brown <dbrown@ucar.edu> wrote:

> Hi Leif,
>
> It can be difficult trying to interpret all the various ways of specifying
> pre-projected grids to be able to plot them in NCL, and this one took
> awhile. But there is always something to learn.
> In this case it is that ESPG (from wikipedia -- the European Petroleum
> Survey Group), which is referenced multiple times in the "spatial_ref"
> attribute has a useful web site that allows you to search for particular
> grids. And that this grid is known on their site as ESPG 3412. Searching on
> their site for 3412 (eventually) leads to this page:
>
> http://spatialreference.org/ref/epsg/3412/
>
> There among other things, you can see the parameters to provide to the
> proj4 map projection tool. Using these parameters I believe I have
> produced an accurate plot of the data -- both
> using the native projection, and using a complete set of lat/lon
> coordinates that I generated, a plot in another projection of the same
> data. I am attaching some png files that demonstrate, along
> with the NCL script I used to produce the plots. I also produced a
> modified version of your data file that contains the lat/lon coordinate
> arrays. This is too big to attach, but you can get it
> from ftp.cgd.ucar.edu as *incoming/leif_ice_llcoords.nc*. I will also
> include the script I used to add the coordinates to the original NetCDF
> file. This depended on creating an ascii file with the
> lat/lon coordinates using 'invproj' first. There is a bit more explanation
> in the comments of the script. Hope this helps.
> -dave
>
>
>
>
>
>
>
>
> On Feb 18, 2013, at 12:17 PM, Leif Truelsen wrote:
>
> Hello,
>
> I have data that in a polar stereographic projection. However, there are
> no associated latitude and longitude variables, so when I attempt to plot
> with ncl I get the errors about not having valid coordinate arrays (e.g.,
> check_for_y_lat_coord: Warning: Data either does not contain a valid
> latitude coordinate array or doesn't contain one at all.). I think the
> most clear solution would be if I created latitude and longitude variables,
> no? Can anyone help direct me how to do this, or of the way to
> plot/manipulate the data as they are?
>
> Best regards,
> Leif
>
>
> Here's my data:
>
> dimensions:
> x = 1779;
> y = 1854;
> time = 10;
> variables:
> char polar_stereographic;
> :grid_mapping_name = "polar_stereographic";
> :straight_vertical_longitude_from_pole = 0.0; // double
> :false_easting = 0.0; // double
> :false_northing = 0.0; // double
> :latitude_of_projection_origin = -90.0; // double
> :standard_parallel = -70.0; // double
> :longitude_of_prime_meridian = 0.0; // double
> :semi_major_axis = 6378273.0; // double
> :inverse_flattening = 298.279411123061; // double
> :spatial_ref = "PROJCS[\"NSIDC Sea Ice Polar Stereographic South\",GEOGCS[\"Unspecified datum based upon the Hughes 1980 ellipsoid\",DATUM[\"Not_specified_based_on_Hughes_1980_ellipsoid\",SPHEROID[\"Hughes 1980\",6378273,298.279411123061,AUTHORITY[\"EPSG\",\"7058\"]],AUTHORITY[\"EPSG\",\"6054\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"4054\"]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",-70],PARAMETER[\"central_meridian\",0],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"X\",EAST],AXIS[\"Y\",NORTH],AUTHORITY[\"EPSG\",\"3412\"]]";
> :GeoTransform = "-3966250 4450 0 4119400 0 -4450 ";
> :_CoordinateTransformType = "Projection";
> :_CoordinateAxisTypes = "GeoX GeoY";
> float ice(time=10, y=1854, x=1779);
> :_FillValue = -3.4028235E38f; // float
> :grid_mapping = "polar_stereographic";
> double time(time=10);
> :calendar = "365_day";
> :standard_name = "time";
> :units = "hours since 2001-00-00 00:00:00";
> double x(x=1779);
> :standard_name = "projection_x_coordinate";
> :long_name = "x coordinate of projection";
> :units = "m";
> :_CoordinateAxisType = "GeoX";
> double y(y=1854);
> :standard_name = "projection_y_coordinate";
> :long_name = "y coordinate of projection";
> :units = "m";
> :_CoordinateAxisType = "GeoY";
>
> _______________________________________________
> 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 Feb 20 13:40:45 2013