[PROJ] GOES Geostationary Projection

classic Classic list List threaded Threaded
6 messages Options
Reply | Threaded
Open this post in threaded view
|

[PROJ] GOES Geostationary Projection

Kurt Schwehr-2

Hi all,

I'm trying to figure out how to convert GOES 16 & 17 weather data into a form readable by GeoTools.  It looks like I'm possibly loosing the sweep information.  Take for example this set of 16 bands:

gsutil cp gs://gcp-public-data-goes-17/ABI-L1b-RadC/2019/151/21/OR_ABI-L1b-RadC-*_G17_s20191512131197_*.nc .

if I do:

gdalwarp --version
GDAL 2.2.4, released 2018/03/19
# e.g. for all C1..C16
gdalwarp -t_srs EPSG:4326 NETCDF:OR_ABI-L1b-RadC-M6C16_G17_s20191512131197_e20191512133581_c20191512134014.nc:Rad OR_ABI-L1b-RadC-M6C16_G17_s20191512131197_e20191512133581_c20191512134014_4326.tif

The imagery once tiled by and displayed by Earth Engine, it shows up in the correct location.

But I would rather have the data stay in it's native projection and only be transformed when used to which ever projection is needed.  (this is typical of Earth Engine)

gdal_translate NETCDF:OR_ABI-L1b-RadC-M6C16_G17_s20191512131197_e20191512133581_c20191512134014.nc:Rad OR_ABI-L1b-RadC-M6C16_G17_s20191512131197_e20191512133581_c20191512134014.tif

Gives the imagery not quite in the correct location.   I briefly looked at gdal_translate from proj at head and gdal at head, but our tiler is an older version of gdal that can read the post "RFC 73: Integration of PROJ6 for WKT2" info.

Any one have suggestions on how to build a geotiff using gdal 2.2.4 from the netcdf that keeps all the necessary projection info in a form that GeoTools 19 can understand?

I've attached a screenshot from the Gulf of California where an island is clearly offset between the renderings of the two assets (transformed from the gdal_translate geotiff and the gdalwarp to EPSG:4326 geotiff)

Thanks!
-kurt
Engineer on Google Earth Engine

Earth Engine screen shots:
https://photos.app.goo.gl/UAGccF2rjPDfenS27

diff -u netcdf.info tif.info
--- netcdf.info 2019-06-26 10:54:22.512939636 -0700
+++ tif.info 2019-06-26 10:54:39.596731603 -0700
@@ -1,22 +1,23 @@
-Driver: netCDF/Network Common Data Format
-Files: OR_ABI-L1b-RadC-M6C02_G17_s20191512131197_e20191512133570_c20191512133595.nc
+Driver: GTiff/GeoTIFF
+Files: OR_ABI-L1b-RadC-M6C02_G17_s20191512131197_e20191512133570_c20191512133595.tif
 Size is 10000, 6000
 Coordinate System is:
-PROJCS["unnamed",
-    GEOGCS["unknown",
-        DATUM["unknown",
+PROJCS["Geostationary_Satellite",
+    GEOGCS["GCS_unknown",
+        DATUM["D_unknown",
             SPHEROID["Spheroid",6378137,298.2572221]],
         PRIMEM["Greenwich",0],
-        UNIT["degree",0.0174532925199433]],
+        UNIT["Degree",0.017453292519943295]],
     PROJECTION["Geostationary_Satellite"],
     PARAMETER["central_meridian",-137],
     PARAMETER["satellite_height",35786023],
     PARAMETER["false_easting",0],
     PARAMETER["false_northing",0],
-    EXTENSION["PROJ4","+proj=geos +lon_0=-137 +h=35786023 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs  +sweep=x"]]
+    UNIT["Meter",1]]
 Origin = (-2505021.493779813405126,1583173.639055640902370)
 Pixel Size = (501.004328871885207,501.004328871885264)
 Metadata:
+  AREA_OR_POINT=Area
   goes_imager_projection#grid_mapping_name=geostationary
   goes_imager_projection#inverse_flattening=298.2572221
   goes_imager_projection#latitude_of_projection_origin=0
@@ -81,23 +82,16 @@
   y#scale_factor=-1.4e-05
   y#standard_name=projection_y_coordinate
   y#units=rad
-Geolocation:
-  LINE_OFFSET=0
-  LINE_STEP=1
-  PIXEL_OFFSET=0
-  PIXEL_STEP=1
-  SRS=GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
-  X_BAND=1
-  X_DATASET=NETCDF:"OR_ABI-L1b-RadC-M6C02_G17_s20191512131197_e20191512133570_c20191512133595.nc":x
-  Y_BAND=1
-  Y_DATASET=NETCDF:"OR_ABI-L1b-RadC-M6C02_G17_s20191512131197_e20191512133570_c20191512133595.nc":y
+Image Structure Metadata:
+  COMPRESSION=DEFLATE
+  INTERLEAVE=BAND
 Corner Coordinates:
-Upper Left  (-2505021.494, 1583173.639) (161d34'43.88"W, 14d47'44.39"N)
-Lower Left  (-2505021.494, 4589199.612) (175d33' 0.44"E, 53d31'38.26"N)
-Upper Right ( 2505021.795, 1583173.639) (112d25'16.11"W, 14d47'44.39"N)
-Lower Right ( 2505021.795, 4589199.612) ( 89d33' 0.41"W, 53d31'38.26"N)
+Upper Left  (-2505021.494, 1583173.639) (161d33'27.65"W, 14d49'58.11"N)
+Lower Left  (-2505021.494, 4589199.612) (175d46'57.22"E, 53d43' 3.62"N)
+Upper Right ( 2505021.795, 1583173.639) (112d26'32.33"W, 14d49'58.11"N)
+Lower Right ( 2505021.795, 4589199.612) ( 89d46'57.19"W, 53d43' 3.62"N)
 Center      (       0.151, 3086186.626) (136d59'59.99"W, 29d58' 0.14"N)
-Band 1 Block=10000x1 Type=Int16, ColorInterp=Undefined
+Band 1 Block=10000x1 Type=Int16, ColorInterp=Gray
   NoData Value=4095
   Unit Type: W m-2 sr-1 um-1
   Offset: -20.2899112701416,   Scale:0.158592373132706

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: GOES Geostationary Projection

Even Rouault-2
> Gives the imagery not quite in the correct location.   I briefly looked at
> gdal_translate from proj at head and gdal at head, but our tiler is an
> older version of gdal that can read the post "RFC 73: Integration of PROJ6
> for WKT2" info.

PROJ6 support in GDAL won't help in any way. This is a pure GeoTIFF issue,
since GeoTIFF doesn't support this projection method, and I'm afraid there is
no clean resolution. Recent enough version of GDALs will incorporate a hacky
WKT in a PCSCitationGeoKey with a EXTENSION["PROJ4", ...] as the example you
shown. But I would be very surprised that GeoTools would understand that !

Actually you might add this use case to this ticket
https://github.com/opengeospatial/geotiff/issues/56 :-) if you want a
ultimately clean solution for that. But even using latest WKT version wouldn't
likely help, since as far as I remember, there is no standardized way in the
EPSG dataset to represent the geostationary projection and the sweep
parameter, so this would involve lobying at that level too.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: GOES Geostationary Projection

Kurt Schwehr-2
Thanks Even!  

Since the end target isn't geotiff and it does have WKT, maybe I can work with the GeoTools folks to get a solution if I cut geotiff out of the middle.

On Wed, Jun 26, 2019 at 3:15 PM Even Rouault <[hidden email]> wrote:
> Gives the imagery not quite in the correct location.   I briefly looked at
> gdal_translate from proj at head and gdal at head, but our tiler is an
> older version of gdal that can read the post "RFC 73: Integration of PROJ6
> for WKT2" info.

PROJ6 support in GDAL won't help in any way. This is a pure GeoTIFF issue,
since GeoTIFF doesn't support this projection method, and I'm afraid there is
no clean resolution. Recent enough version of GDALs will incorporate a hacky
WKT in a PCSCitationGeoKey with a EXTENSION["PROJ4", ...] as the example you
shown. But I would be very surprised that GeoTools would understand that !

Actually you might add this use case to this ticket
https://github.com/opengeospatial/geotiff/issues/56 :-) if you want a
ultimately clean solution for that. But even using latest WKT version wouldn't
likely help, since as far as I remember, there is no standardized way in the
EPSG dataset to represent the geostationary projection and the sweep
parameter, so this would involve lobying at that level too.

Even

--
Spatialys - Geospatial professional services
http://www.spatialys.com


--

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: GOES Geostationary Projection

Even Rouault-2
On jeudi 27 juin 2019 08:43:06 CEST Kurt Schwehr wrote:
> Thanks Even!
>
> Since the end target isn't geotiff and it does have WKT

Which WKT are you talking about ?
I consider the WKT1 with PROJ strings you'll find reported by GDAL on the
netCDF or in GeoTIFF to be hack, but I don't know a standardized form of WKT
for that.

With GDAL 3 / PROJ 6, the WKT2 is probably a bit better (see just below). But
the method name and parameter names are still my own invention (hopefully
reasonable enough)

$ gdalsrsinfo "+proj=geos +sweep=x +h=35786023"

PROJ.4 : +proj=geos +sweep=x +lon_0=0 +h=35786023 +x_0=0 +y_0=0 +datum=WGS84
+units=m +no_defs

OGC WKT2:2018 :
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Geostationary Satellite (Sweep X)"],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Satellite Height",35786023,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]


--
Spatialys - Geospatial professional services
http://www.spatialys.com
_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: GOES Geostationary Projection

Kurt Schwehr-2
The GDAL in my tool chain is circa Oct-2018, so the hack.  But since I'm writing a custom converter, I should be able to force the WKT to be what ever GeoTools wants for GEOS 16 & 17.  I'll start up a discussion on the GeoTools mailing list and link back to here.

On Thu, Jun 27, 2019 at 9:32 AM Even Rouault <[hidden email]> wrote:
On jeudi 27 juin 2019 08:43:06 CEST Kurt Schwehr wrote:
> Thanks Even!
>
> Since the end target isn't geotiff and it does have WKT

Which WKT are you talking about ?
I consider the WKT1 with PROJ strings you'll find reported by GDAL on the
netCDF or in GeoTIFF to be hack, but I don't know a standardized form of WKT
for that.

With GDAL 3 / PROJ 6, the WKT2 is probably a bit better (see just below). But
the method name and parameter names are still my own invention (hopefully
reasonable enough)

$ gdalsrsinfo "+proj=geos +sweep=x +h=35786023"

PROJ.4 : +proj=geos +sweep=x +lon_0=0 +h=35786023 +x_0=0 +y_0=0 +datum=WGS84
+units=m +no_defs

OGC WKT2:2018 :
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Geostationary Satellite (Sweep X)"],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Satellite Height",35786023,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]


--
Spatialys - Geospatial professional services
http://www.spatialys.com


--

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj
Reply | Threaded
Open this post in threaded view
|

Re: GOES Geostationary Projection

Kurt Schwehr-2

On Thu, Jun 27, 2019 at 10:08 AM Kurt Schwehr <[hidden email]> wrote:
The GDAL in my tool chain is circa Oct-2018, so the hack.  But since I'm writing a custom converter, I should be able to force the WKT to be what ever GeoTools wants for GEOS 16 & 17.  I'll start up a discussion on the GeoTools mailing list and link back to here.

On Thu, Jun 27, 2019 at 9:32 AM Even Rouault <[hidden email]> wrote:
On jeudi 27 juin 2019 08:43:06 CEST Kurt Schwehr wrote:
> Thanks Even!
>
> Since the end target isn't geotiff and it does have WKT

Which WKT are you talking about ?
I consider the WKT1 with PROJ strings you'll find reported by GDAL on the
netCDF or in GeoTIFF to be hack, but I don't know a standardized form of WKT
for that.

With GDAL 3 / PROJ 6, the WKT2 is probably a bit better (see just below). But
the method name and parameter names are still my own invention (hopefully
reasonable enough)

$ gdalsrsinfo "+proj=geos +sweep=x +h=35786023"

PROJ.4 : +proj=geos +sweep=x +lon_0=0 +h=35786023 +x_0=0 +y_0=0 +datum=WGS84
+units=m +no_defs

OGC WKT2:2018 :
PROJCRS["unknown",
    BASEGEOGCRS["unknown",
        DATUM["World Geodetic System 1984",
            ELLIPSOID["WGS 84",6378137,298.257223563,
                LENGTHUNIT["metre",1]],
            ID["EPSG",6326]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8901]]],
    CONVERSION["unknown",
        METHOD["Geostationary Satellite (Sweep X)"],
        PARAMETER["Longitude of natural origin",0,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8802]],
        PARAMETER["Satellite Height",35786023,
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        PARAMETER["False easting",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8806]],
        PARAMETER["False northing",0,
            LENGTHUNIT["metre",1],
            ID["EPSG",8807]]],
    CS[Cartesian,2],
        AXIS["(E)",east,
            ORDER[1],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]],
        AXIS["(N)",north,
            ORDER[2],
            LENGTHUNIT["metre",1,
                ID["EPSG",9001]]]]


--
Spatialys - Geospatial professional services
http://www.spatialys.com


--


--

_______________________________________________
PROJ mailing list
[hidden email]
https://lists.osgeo.org/mailman/listinfo/proj