Add Lambert Conformal Conic single standard parallel support to netCDF driver

[kyle]

NetCDF driver does not support a single standard parallel for lambert conformal conic, and it sets the two standard parallel version to the default 0, even if a single parallel exists. I suggest adding support Lambert_Conformal:standard_parallel tag for a single standard parallel and possibly check and see if Lambert_Conformal:standard_parallel_1 == Lambert_Conformal:standard_parallel_2 as well.

Attached is a netcdf file that exhibits the behavior

[kyle]

netcdf dataset

[kyle]

According to the CF-1.0 convention I found a few problems in the driver. I don't know if the netcdf driver should follow CF strictly or not though. There were problems retrieving the values associated with the following netcdf tags:

standard_parallel_1 standard_parallel_2 central_meridian

None of these are in CF-1.0

In the CF-1.0 convention there is no distinction between the tags for standard parallels. The tag is: standard_parallel and it may occur either once or twice.

The CF-1.0 tag name for the central meridian is: longitude_of_central_meridian

There are a couple more issues, but I don't think they fall under the scope of this ticket, as one is a units issue and the other is a projection issue with the datum (the datum for most projections defaults to WGS-84 even if the definition of the sphere is incuded in the netcdf file)

I am slowly working on a patch, although my solution will surely not be as elegant as others.

link to cf-1.0 convention grid_mapping tags:  http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.4/cf-conventions.html#id3000051

Thanks, kss

[kyle]

My reference is CF-1.4, not CF-1.0, sorry.

[kyle]

lcc 1sp patch

[kyle]

Adding a patch that specifically supports standard_parallel tag, it is only called in the lcc section of SetProjection?(), it still uses the SetLCC() method of OGRSpatialReference and sets both parallels equal to each other. This is okay as long as the latitude of origin is set as well.

[warmerdam]

Kyle,

I'm getting:

warmerda@gdal64[75]% patch --dry-run < ~/Desktop/add_lcc1sp_netcdfdataset_kyle.patch
patching file netcdfdataset.h
Reversed (or previously applied) patch detected!  Assume -R? [n] n
Apply anyway? [n] y
Hunk #1 FAILED at 58.
Hunk #2 succeeded at 162 (offset 10 lines).
1 out of 2 hunks FAILED -- saving rejects to file netcdfdataset.h.rej
patching file netcdfdataset.cpp
Hunk #4 succeeded at 935 with fuzz 1.

Perhaps this overlaps with the patch for #3325?

[kyle]

Yes, it does, I am working on a fix.

[kyle]

r18896 checks cf and gdal standard_parallel tags, but not cf central_meridian tags

[kyle]

added lcc1sp support and updated in trunk (r18896)

[etiennesky]

Sorry to re-open this old ticket, but I have just noticed that the driver generates a "lambert_conformal_conic1" grid mapping name, which is not a valid CF grid mapping type.

I would suggest to set is as a "lambert_conformal_conic" grid mapping, type with only one standard_parallel tag. This would make it compatible with other software, and also for input into GDAL, since reading has been fixed.

I am working on a patch for fixing incorrect standard_parallel_1 and standard_parallel_2 attributes (issue #2893), and would like to include this as well.

Sounds reasonable? Perhaps this doesn't belong in this ticket, but it is highly related.

Found this issue with the cf_lcc1sp.nc file of the autotests. Same procedure with double standard parallel (cf_lcc2sp.nc) generates valid CF-1.x.

gdal_translate -of netcdf cf_lcc1sp.nc cf_lcc1sp-mod1.nc

ncdump -h cf_lcc1sp-mod1.nc

netcdf cf_lcc1sp-mod1 { dimensions:

x = 11 ; y = 14 ;

variables:

char lambert_conformal_conic1 ;

...

lambert_conformal_conic1:grid_mapping_name = "lambert_conformal_conic1" ; lambert_conformal_conic1:latitude_of_projection_origin = 25. ;

... }

cfchecks -v 1.0 cf_lcc2sp-mod1.nc

... ERROR (5.6): Invalid grid_mapping_name: lambert_conformal_conic1 ...

[kyle]

Etienne, I agree this needs to be handled, possible a new ticket for writing correct cf convention projections. That would limit the scope a little and make it easier to manage. I didn't do any checks on writing when I implemented the read part. I don't know if there are any other projections that are written incorrectly, but that is easy to check. I will also make some comments on the wiki page when I get a chance.

kss

[etourigny]

Can I safely use oSRS.SetLCC() when there are 1 or 2 standard parallels? Internally, does SetLCC1SP() amount to the same as SetLCC(StdP1, StdP1) - in the case of one parallel StdP1?

This would be a very simple fix, but I want it to be correct.

In terms of CF, there is no difference between the two, and they are all called lambert_conformal_conic.

[rouault]

I'm definitely not the reference on that topic (Frank would be of better advice). Perhaps  http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html and  http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_2sp.html ( reached from the home page of  http://trac.osgeo.org/geotiff/ : Catalog of projections link ) might give some clues. Grep'ping through GDAL source tree could also bring some light as I think there are other drivers where the distinction between 1SP and 2SP is not clearly made.

[etourigny]

Fixed the export in r23265 so that standard_parallel_1 and standard_parallel_2 are exported to standard_parallel with 2 values.

However, LCC1SP is still not fully resolved (as LCC does not have scale_factor in CF-1), see wiki:NetCDF_ProjectionTestingStatus.

[etourigny]

in trunk (r23579): fix LCC-1SP import and export(bug #3324) - added computation of scale factor if it is missing

I have come up with a reasonable solution to manage the inconsistencies between WKT and CF for the 1SP variant. More info at  http://mailman.cgd.ucar.edu/pipermail/cf-metadata/2011/022650.html ,  http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html and  http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.5/cf-conventions.html#appendix-grid-mappings

CF uses standard_parallel1 and latitude_of_projection_origin but no scale factor. If they are equal (the common case) then scale factor is 1.0 , else one can obtain it with Snyder eq. 15-4. The reverse computation has no solution (that I know of), so the standard_parallel1 parameter is kept and exported for later use (as it is CF). In the case that scale factor !=1 and standard_parallel1 is not available (standard OGC WKT), then an error is thrown and using the 2SP variant is suggested.

I have tested the patch with the file provided here and it works well, although I did not test the implementation of eq. 15-4 fully. I would need more test cases with actual data (and known scale factor).

Ideally, the 1SP variant with scale factor != 1 should be transformed to a 2SP form, but I don't know the formulas for that.

Kyle: If you follow this, can you give your opinion and perhaps test it?

[etourigny]

in trunk (r23580): fix LCC-1SP export(bug #3324) - if std_parallel is missing and scale factor=1, copy it from latitude_of_origin

This is the converse of setting scale_factor=1 when latitude_of_origin==std_parallel on import.

[kyle]

Etienne, I am getting peculiar results with the new lcc 1sp in 1.9: 1.8.1:

Corner Coordinates:
Upper Left  (-1696476.103, 2404012.657) (115d16'30.26"W, 45d 2' 4.70"N)
Lower Left  (-1696476.103, 2099237.657) (114d44'54.00"W, 42d28'34.77"N)
Upper Right (-1452656.103, 2404012.657) (112d23'44.19"W, 45d19' 6.52"N)
Lower Right (-1452656.103, 2099237.657) (111d56'31.04"W, 42d45'27.03"N)
Center      (-1574566.103, 2251625.157) (113d35'20.71"W, 43d54'27.73"N)

1.9dev:

Corner Coordinates:
Upper Left  (-1696476.103, 2404012.657) ( 77d14'16.21"W, 89d59'58.03"S)
Lower Left  (-1696476.103, 2099237.657) ( 68d20'10.10"W, 89d59'57.58"S)
Upper Right (-1452656.103, 2404012.657) ( 86d52'59.27"W, 89d59'57.80"S)
Lower Right (-1452656.103, 2099237.657) ( 78d26' 5.12"W, 89d59'57.24"S)
Center      (-1574566.103, 2251625.157) ( 77d47'37.31"W, 89d59'57.68"S)

These files *seemed* to be correct when I implemented lcc the first time. I will try to look into some lcc documentation.

file attached.

[kyle]

Also, I cannot warp this file any more: 1.8.1:

[kyle@localhost Desktop]$ gdalwarp -t_srs EPSG:4326 NETCDF:20110401T1800.nc:Total_cloud_cover 181.tif
Creating output file that is 26P x 23L.
Processing input file NETCDF:20110401T1800.nc:Total_cloud_cover.
Using internal nodata values (eg. nan) for image NETCDF:20110401T1800.nc:Total_cloud_cover.
0...10...20...30...40...50...60...70...80...90...100 - done.

1.9dev:

[kyle@localhost Desktop]$ gdalwarp -t_srs EPSG:4326 NETCDF:20110401T1800.nc:Total_cloud_cover 19dev.tif
Creating output file that is 496P x -2147483648L.
ERROR 1: Attempt to create 496x-2147483648 dataset is illegal,sizes must be larger than zero.

[kyle]

The offending code is here:

/* store dfStdP1 so we can output it to CF later */
                        oSRS.SetNormProjParm( SRS_PP_STANDARD_PARALLEL_1, dfStdP1 );
                    }

If a standard parallel is stored, it is used in projections. Looking at ogr_spatialref.h, LCC1SP doesn't have any standard parallels. I may just have a bad file. Does CF have a tag for latitude of origin?

[etourigny]

I have been trying (with some prior help from Pat Sunter) to make the CF<->WKT translation work well with LCC1SP, the problem is that the CF definition does not match the OGR WKT definition.

Please see http://trac.osgeo.org/gdal/wiki/NetCDF_ProjectionTestingStatus  http://cf-pcmdi.llnl.gov/documents/cf-conventions/1.5/cf-conventions.html#appendix-grid-mappings  http://www.remotesensing.org/geotiff/proj_list/lambert_conic_conformal_1sp.html

and read this thread (there are 3 projections that are discussed)

 http://mailman.cgd.ucar.edu/pipermail/cf-metadata/2011/027605.html

as you can see from the CF-1.5 page, there IS a latitude_of_projection_origin but no scale_factor. In the case that k0=1, the translation WKT<->CF is done properly

Lambert conformal
grid_mapping_name = lambert_conformal_conic
        
Map parameters:
standard_parallel - There may be 1 or 2 values.
longitude_of_central_meridian
latitude_of_projection_origin
false_easting
false_northing

Pat Sunter's solution of some time ago (to keep standard parallel 1) seemed reasonable, as CF uses that parameter. I didn't realise it was used in the projection though, as it is a LCC-1SP in OGR... Can you confirm that if you remove that one line the reprojection and warping works as before?

The fix in r23579 and r23580 was to relate standard_parallel and scale_factor. I guess I should remove standard_parallel entirely now from the WKT in OGR. I'll try that.

Thanks for your help

[etourigny]

r23611 : fix handling of UNITS import and export (bugs #4402 and #3324)

The bug was in due to the double scaling problem described in #4402 . The standard_parallel parameter does not affect the transformation, I checked the proj4 string is the same. In fact the difference was in the +units value (m vs. km). This does not surprise me, because LCC-1SP does not use std_parallel for LCC-1SP. I will keep that parameter for now, because it is needed to preserve CF tags. See the code marked "special case for LCC-1SP" in netcdfdataset.cpp for more details.

Please confirm this works for you also. Thanks again

[kyle]

Looks good to me, thanks etienne.

[etourigny]

see #2893 for more details