[R-sig-Geo] NASA's Black Marble monthly data: Reprojection isn't accurate

[Barry Rowlingson]

My gdal 3.4.1 doesn't get the extent.

I answered this (badly) on gis.stackexchange - its a 10degree tile where
the origin is related to the h and v elements of the path (and also stored
as attributes in the netcdf structure).

I think I messed up the vertical offset, which was because I only bothered
making it work on one tile. The tile position is going to be A+Bh, C+Dv
where A,B,C,D are constants, and B and D are either plus or minus 10
depending on if the origin is top or bottom or left or right, and A and C
are the origins of the 0-th tiles. Once you have the origin, add 10 to get
the extent. This gives the correct resolution (15 minutes of arc).

Maybe tomorrow I'll get a few more tiles and write up my answer in an edit
on stackex, and check all the extents etc are correct, and put my code in a
function. I didn't have time to do all that previously but hoped there was
enough there for understanding of the problem and anyone with a bit of
maths could fix it up.

Barry


On Sat, Sep 9, 2023 at 9:28 PM Michael Sumner <mdsumner using gmail.com> wrote:

> The extent of this one is
>
> 70, 80, 10, 20
>
> Later versions of GDAL determine this automatically (I'm not sure when)
>
>  rast("VNP46A3.A2018091.h25v07.001.2021125122857.h5")
> class       : SpatRaster
> dimensions  : 2400, 2400, 26  (nrow, ncol, nlyr)
> resolution  : 0.004166667, 0.004166667  (x, y)
> extent      : 70, 80, 10, 20  (xmin, xmax, ymin, ymax)
> coord. ref. : lon/lat Unknown datum based upon the GRS 1980 Authalic Sphere
> ellipsoid
> sources     :
>
> VNP46A3.A2018091.h25v07.001.2021125122857.h5://AllAngle_Composite_Snow_Covered
>
>
> VNP46A3.A2018091.h25v07.001.2021125122857.h5://AllAngle_Composite_Snow_Covered_Num
>
>
> VNP46A3.A2018091.h25v07.001.2021125122857.h5://AllAngle_Composite_Snow_Covered_Quality
>               ... and 23 more source(s)
> varnames    : AllAngle_Composite_Snow_Covered
>               AllAngle_Composite_Snow_Covered_Num
>               AllAngle_Composite_Snow_Covered_Quality
>               ...
>
> Be very careful with canned steps to "fix" georeferencing, make sure it
> needs fixing and check that it's right. Ultimately the best way to "fix" it
> is to go to the source of the code that interfaces the data, which here is
> GDAL and report there - but clearly that's been updated from whatever
> version was in use on stackoverflow.
>
> Cheers, Mike
> .
>
>
>
>
>
> On Sun, Sep 10, 2023 at 1:17 AM Nikolaos Tziokas <nikos.tziokas using gmail.com>
> wrote:
>
> > I downloaded NASA's Black Marble monthly nighttime light NTL data,
> VNP46A3
> > <
> >
> https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A3/
> > >.
> > In a previous question
> > <
> >
> https://gis.stackexchange.com/questions/466571/extent-not-found-on-nasas-black-marble-monthly-images-how-to-set-it/466574?noredirect=1#comment761916_466574
> > >
> > of mine the reprojection worked perfectly but now it seems that it
> doesn't.
> > For example, I wanted to download NTL data for the city of Mumbai, India.
> > After reprojecting the NTL (product 5 (All_Angles_Snow_Free) from the
> .h5)
> > the result is attached here
> > <
> https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7
> > >.
> >
> > At the bottom if the image is a shp of Mumbai (downloaded from GADM) and
> > the red circle in the top indicates where Mumbai is in the NTL image.
> > Clearly something's not right.
> >
> > I downloaded the image from here
> > <
> >
> https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/VNP46A3/
> > >
> > (LAADS-DAAC, Level-1 and Atmosphere Archive & Distribution System
> > Distributed Active Archive Center). The code I used to extract the NTL
> > radiance image is:
> >
> > library(terra)
> >
> > wd <- "path/"
> >
> > r <- rast(paste0(wd, "VNP46A3.A2018091.h25v07.001.2021125122857.h5"))
> > crs(r) <- "epsg:4326"
> >
> > 2400*(15/(60*60))
> >
> > h = 25
> > v = 7
> >
> > ext(r) = c(-180+h*10,-180+(h+1)*10, (v-2)*10,(v-1)*10)
> >
> > ntl <- r[[5]]
> >
> > writeRaster(ntl, paste0(wd, "ntl.tif"), overwrite = TRUE)
> >
> > Why the code worked perfectly in the previous question and now it
> doesn't?
> > From here
> > <
> https://drive.google.com/drive/folders/1V115zpdU2-5fXssI6iWv_F6aNu4E5qA7>
> > you can download the .h5 image if you don't want to use NASA's website. I
> > am using R 4.3.1 and RStudio 2023.06.2+561.
> >
> > --
> > Tziokas Nikolaos
> > Cartographer
> >
> > Tel:(+44)07561120302
> > LinkedIn <http://linkedin.com/in/nikolaos-tziokas-896081130>
> >
> >         [[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-Geo mailing list
> > R-sig-Geo using r-project.org
> > https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >
>
>
> --
> Michael Sumner
> Software and Database Engineer
> Australian Antarctic Division
> Hobart, Australia
> e-mail: mdsumner using gmail.com
>
>         [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo using r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]