[ncl-talk] Projection lat/lon (follow-up to earlier GeoTIFF to netcdf email)
Rick Brownrigg
brownrig at ucar.edu
Thu Mar 8 13:51:38 MST 2018
Hi Michael,
Here are two solutions to converting your geotiffs into geographic
coordinates. The first one uses GDAL, and is a kind of hackish/one-shot
approach that involves hand-jamming some number in; its virtue is that its
easy, but would not scale if you had to do convert many files. The second
approach uses NCL to do the conversion, and tries to take a programmatic
approach to extracting the projection parameters from the file
i) I think you mentioned using gdal_translate to convert your tiff into
NetCDF. The GDAL utility "gdalwarp" will do that conversion and perform the
reprojection. I used:
gdalwarp -s_srs "+proj=lcc +lon_0=-100 +lat_0=42.5 +lat_1=25 +lat_2=60
+a=6378137 +b=6356752" -t_srs "+proj=longlat" -of netCDF
slope_ncols_nrows_rescale_flip_shift.tif
slope_ncols_nrows_rescale_flip_shift.nc
The s_srs and t_srs parameters are source/target projection parameters.
They are given whats referred to as a "proj4 string", which is the manner
in which the Proj4 cartographic projection library expects projection
parameters. The numerical values in the s_srs where extracted by looking at
the metadata of the original geotiff. Note that the "-of netCDF" specifies
the output file format.
Attached is an NCL script that simply reads the reprojected file and plots
it (no attempt at esthetics here -- I just want to know of the
georeferencing was correct!). (plotslope.ncl and plotslope.png).
ii) The second approach still uses gdal_translate to convert geotiff into
NetCDF. It then utilizes NCL's (undocumented) interface to the Proj4
library. It first programmatically extracts metadata from the NetCDF to
build up proj4 strings and then invokes transform_coordinate() to convert
the lambert conformal coordinates into lat/lon. It also produces a plot
just to verify that the georeferencing is correct. (xy2LL.ncl and
xy2LL.png). Note that this script was rather memory intensive, taking more
memory resource than my local 16GB machine could offer.
So, I hope that helps. Let me know if there are questions about any of this.
Rick
On Fri, Mar 2, 2018 at 8:13 AM, Michael Notaro <mnotaro at wisc.edu> wrote:
> That's fine, I can wait. I appreciate your help.
>
> Michael
>
>
> Michael Notaro
> Associate Director
> Nelson Institute Center for Climatic Research
> University of Wisconsin-Madison
> Phone: (608) 261-1503
> Email: mnotaro at wisc.edu
>
>
> ------------------------------
> *From:* Rick Brownrigg <brownrig at ucar.edu>
> *Sent:* Friday, March 2, 2018 8:31 AM
>
> *To:* Michael Notaro
> *Subject:* Re: [ncl-talk] Projection lat/lon (follow-up to earlier
> GeoTIFF to netcdf email)
>
> Hi Michael,
>
> I was able to get those files. Unfortunately, I won't be able to do much
> with them until Monday evening at the earliest, as I'm on travel starting
> today. I don't know you timeline, but feel free to re-post to ncl-talk if
> things are more urgent for you.
>
> Rick
>
>
> On Thu, Mar 1, 2018 at 1:59 PM, Michael Notaro <mnotaro at wisc.edu> wrote:
>
> Thanks. I finished ftp'ing the files. The filenames are:
>
>
> slope_daymet.nc
>
>
> slope_ncols_nrows_rescale_flip_shift.tfw
>
> slope_ncols_nrows_rescale_flip_shift.tif
>
> slope_ncols_nrows_rescale_flip_shift.tif.aux.xml
>
> slope_ncols_nrows_rescale_flip_shift.tif.ovr
>
> slope_ncols_nrows_rescale_flip_shift.tif.xml
>
> aspect_daymet.nc
>
>
> aspect_ncols_nrows_rescale_flip_shift.tfw
>
> aspect_ncols_nrows_rescale_flip_shift.tif
>
> aspect_ncols_nrows_rescale_flip_shift.tif.aux.xml
>
> aspect_ncols_nrows_rescale_flip_shift.tif.ovr
>
> aspect_ncols_nrows_rescale_flip_shift.tif.xml
>
> Michael
>
>
> Michael Notaro
> Associate Director
> Nelson Institute Center for Climatic Research
> University of Wisconsin-Madison
> Phone: (608) 261-1503
> Email: mnotaro at wisc.edu
>
>
> ------------------------------
> *From:* Rick Brownrigg <brownrig at ucar.edu>
> *Sent:* Thursday, March 1, 2018 2:52 PM
> *To:* Michael Notaro
>
> *Subject:* Re: [ncl-talk] Projection lat/lon (follow-up to earlier
> GeoTIFF to netcdf email)
>
> Yes, we an FTP site:
>
> ftp.cgd.ucar.edu
> cd incoming
>
> Please let me know the exact name of the file once posted -- we can't view
> the contents of that directory for "security reasons", although we can
> retreive anything from there if we know the filename.
>
> RB
>
> On Thu, Mar 1, 2018 at 1:48 PM, Michael Notaro <mnotaro at wisc.edu> wrote:
>
> Thanks, Rick. Much appreciated. Do you have an ftp site to
>
> provide the data? Otherwise I can put the files on our own ftp.
>
> Michael
>
>
> Michael Notaro
> Associate Director
> Nelson Institute Center for Climatic Research
> University of Wisconsin-Madison
> Phone: (608) 261-1503
> Email: mnotaro at wisc.edu
>
>
> ------------------------------
> *From:* Rick Brownrigg <brownrig at ucar.edu>
> *Sent:* Thursday, March 1, 2018 2:37 PM
> *To:* Michael Notaro; Ncl-talk
> *Subject:* Re: [ncl-talk] Projection lat/lon (follow-up to earlier
> GeoTIFF to netcdf email)
>
> Hi Michael,
>
> NCL has an unadvertised interface to the Proj4 cartographic library, which
> could be used to convert between the implied meters coordinate system of
> the projection into lat/lons. All the info is there in the metadata needed
> to make this possible.
>
> The function is TransformCoordinate_W(SrcProj:string, DstProj:string,
> x:double, y:double, z:double)
>
> where Src/DstProj are proj4 projection parameters. Coordinates are passed
> in and returned via the x/y/z arrays.
>
> If you are not familiar with Proj4, if you could share this NetCDF file
> with me, I'll take a stab at an implementation showing the steps, and
> repost the solution back to you and as well as the email group. Besides
> that, the original tiff appears to contain multiple subsampled images of
> the original -- I'm curious as to what GDAL did with that.
>
> That all said, I'm traveling the next 4 days, and may not be able to get
> to this right away.
>
> FWIW...
> Rick
>
>
> On Thu, Mar 1, 2018 at 9:29 AM, Michael Notaro <mnotaro at wisc.edu> wrote:
>
> Thanks to Jonathan's suggestion, I used GDAS to convert a GeoTIFF file
> containing slope data to netcdf.
>
> The netcdf file contains Band1(y,x) where y=8732 and x=10405. The
> projection is Lambert_Conformal_Conic_2SP.
>
> See the info below on the projection. How can I obtain the latitude and
> longitude
>
> values for each grid cell?
>
>
> Thanks, Michael
>
>
>
> NCDUMP RESULTS
>
>
> [mmia:~/nps_daymet/slope_ncols_nrows_rescale_flip_shift.tif] notaro%
> ncdump -h slope_daymet.nc
>
> netcdf slope_daymet {
>
> dimensions:
>
> x = 10405 ;
>
> y = 8732 ;
>
> variables:
>
> char lambert_conformal_conic ;
>
> lambert_conformal_conic:grid_mapping_name = "lambert_conformal_conic" ;
>
> lambert_conformal_conic:longitude_of_central_meridian = -100. ;
>
> lambert_conformal_conic:false_easting = 0. ;
>
> lambert_conformal_conic:false_northing = 0. ;
>
> lambert_conformal_conic:latitude_of_projection_origin = 42.5 ;
>
> lambert_conformal_conic:standard_parallel = 25., 60. ;
>
> lambert_conformal_conic:long_name = "CRS definition" ;
>
> lambert_conformal_conic:longitude_of_prime_meridian = 0. ;
>
> lambert_conformal_conic:semi_major_axis = 6378137. ;
>
> lambert_conformal_conic:inverse_flattening = 298.257223563 ;
>
> lambert_conformal_conic:spatial_ref = "PROJCS[\"Lambert_Conformal_Conic_2SP\",GEOGCS[\"WGS
> 84\",DATUM[\"WGS_1984\",SPHEROID[\"WGS 84\",6378137,298.257223563,AUT
> HORITY[\"EPSG\",\"7030\"]],AUTHORITY[\"EPSG\",\"6326\"]],PRI
> MEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433],AUT
> HORITY[\"EPSG\",\"4326\"]],PROJECTION[\"Lambert_Conformal_Co
> nic_2SP\"],PARAMETER[\"standard_parallel_1\",25],PARAMETER[\
> "standard_parallel_2\",60],PARAMETER[\"latitude_of_origin\",
> 42.5],PARAMETER[\"central_meridian\",-100],PARAMETER[\"
> false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"
> metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]" ;
>
> lambert_conformal_conic:GeoTransform = "-5813750 1000 0 4995500 0 -1000 "
> ;
>
> double x(x) ;
>
> x:standard_name = "projection_x_coordinate" ;
>
> x:long_name = "x coordinate of projection" ;
>
> x:units = "m" ;
>
> double y(y) ;
>
> y:standard_name = "projection_y_coordinate" ;
>
> y:long_name = "y coordinate of projection" ;
>
> y:units = "m" ;
>
> float Band1(y, x) ;
>
> Band1:long_name = "GDAL Band Number 1" ;
>
> Band1:_FillValue = -9999.f ;
>
> Band1:SourceBandIndex = 0 ;
>
> Band1:grid_mapping = "lambert_conformal_conic" ;
>
>
> // global attributes:
>
> :GDAL_AREA_OR_POINT = "Area" ;
>
> :GDAL_DataType = "Generic" ;
>
> :Conventions = "CF-1.5" ;
>
> :GDAL = "GDAL 2.2.3, released 2017/11/20" ;
>
> :history = "Thu Mar 01 10:05:20 2018: GDAL CreateCopy( slope_daymet.nc,
> ... )" ;
>
>
> INFO FROM GDAL
>
>
> [mmia:~/nps_daymet/slope_ncols_nrows_rescale_flip_shift.tif] notaro%
> gdalinfo slope_ncols_nrows_rescale_flip_shift.tif
>
> Driver: GTiff/GeoTIFF
>
> Files: slope_ncols_nrows_rescale_flip_shift.tif
>
> slope_ncols_nrows_rescale_flip_shift.tif.ovr
>
> slope_ncols_nrows_rescale_flip_shift.tif.aux.xml
>
> Size is 10405, 8732
>
> Coordinate System is:
>
> PROJCS["Lambert_Conformal_Conic_2SP",
>
> GEOGCS["WGS 84",
>
> DATUM["WGS_1984",
>
> SPHEROID["WGS 84",6378137,298.257223563,
>
> AUTHORITY["EPSG","7030"]],
>
> AUTHORITY["EPSG","6326"]],
>
> PRIMEM["Greenwich",0],
>
> UNIT["degree",0.0174532925199433],
>
> AUTHORITY["EPSG","4326"]],
>
> PROJECTION["Lambert_Conformal_Conic_2SP"],
>
> PARAMETER["standard_parallel_1",25],
>
> PARAMETER["standard_parallel_2",60],
>
> PARAMETER["latitude_of_origin",42.5],
>
> PARAMETER["central_meridian",-100],
>
> PARAMETER["false_easting",0],
>
> PARAMETER["false_northing",0],
>
> UNIT["metre",1,
>
> AUTHORITY["EPSG","9001"]]]
>
> Origin = (-5813750.000000000000000,4995500.000000000000000)
>
> Pixel Size = (1000.000000000000000,-1000.000000000000000)
>
> Metadata:
>
> AREA_OR_POINT=Area
>
> DataType=Generic
>
> Image Structure Metadata:
>
> COMPRESSION=LZW
>
> INTERLEAVE=BAND
>
> Corner Coordinates:
>
> Upper Left (-5813750.000, 4995500.000) (150d36'29.26"E, 47d26' 3.93"N)
>
> Lower Left (-5813750.000,-3736500.000) (142d57' 5.67"W, 2d31'41.12"S)
>
> Upper Right ( 4591250.000, 4995500.000) ( 3d59'18.41"E, 58d17'42.04"N)
>
> Lower Right ( 4591250.000,-3736500.000) ( 64d56'47.60"W, 1d29'22.14"N)
>
> Center ( -611250.000, 629500.000) (108d35'55.06"W, 48d 8'15.05"N)
>
> Band 1 Block=128x128 Type=Float32, ColorInterp=Gray
>
> Min=0.000 Max=0.822
>
> Minimum=0.000, Maximum=0.822, Mean=0.008, StdDev=0.030
>
> NoData Value=-9999
>
> Overviews: 5203x4366, 2602x2183, 1301x1092, 651x546, 326x273, 163x137
>
> Metadata:
>
> SourceBandIndex=0
>
> STATISTICS_COVARIANCES=0.0009288078976163249
>
> STATISTICS_MAXIMUM=0.82241380214691
>
> STATISTICS_MEAN=0.0081764099199458
>
> STATISTICS_MINIMUM=0
>
> STATISTICS_SKIPFACTORX=1
>
> STATISTICS_SKIPFACTORY=1
>
> STATISTICS_STDDEV=0.030476349807947
>
>
> _______________________________________________
> ncl-talk mailing list
> ncl-talk at ucar.edu
> List instructions, subscriber options, unsubscribe:
> http://mailman.ucar.edu/mailman/listinfo/ncl-talk
> ncl-talk Info Page - University Corporation for ...
> <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
> mailman.ucar.edu
> This email list is for NCL users to exchange information, ask questions,
> and report bugs on NCL or other related software if appropriate. It is also
> for the ...
>
> ncl-talk Info Page - University Corporation for ...
> <http://mailman.ucar.edu/mailman/listinfo/ncl-talk>
> mailman.ucar.edu
> Mailing List Manager <http://mailman.ucar.edu/>
> mailman.ucar.edu
> Mailing List Manager. Are you coming here as: Mailing List User. or
> Mailing List Manager/Administrator (requires list pass words) WARNING TO
> LIST ADMINISTRATORS: Do ...
>
> This email list is for NCL users to exchange information, ask questions,
> and report bugs on NCL or other related software if appropriate. It is also
> for the ...
>
>
>
>
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180308/b6c0ad3f/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plotslope.ncl
Type: application/octet-stream
Size: 1038 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180308/b6c0ad3f/attachment-0002.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: xy2LL.ncl
Type: application/octet-stream
Size: 1366 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180308/b6c0ad3f/attachment-0003.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: plotslope.png
Type: image/png
Size: 82977 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180308/b6c0ad3f/attachment-0002.png>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: xy2LL.png
Type: image/png
Size: 108227 bytes
Desc: not available
URL: <http://mailman.ucar.edu/pipermail/ncl-talk/attachments/20180308/b6c0ad3f/attachment-0003.png>
More information about the ncl-talk
mailing list