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

2023-09-10 Thread 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  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.00417, 0.00417  (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 
> 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 
> >
> > [[alternative HTML version deleted]]
> >
> > ___
> > R-sig-Geo mailing list
> > R-sig-Geo@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: mdsum...@gmail.com
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@r-project.org
>

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

2023-09-09 Thread Michael Sumner 
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.00417, 0.00417  (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 
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
>  >.
>
> 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
> 
> 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 
>
> [[alternative HTML version deleted]]
>
> ___
> R-sig-Geo mailing list
> R-sig-Geo@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: mdsum...@gmail.com

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

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

2023-09-09 Thread Nikolaos Tziokas 
I downloaded NASA's Black Marble monthly nighttime light NTL data, VNP46A3
.
In a previous question

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
.

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

(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

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 

[[alternative HTML version deleted]]

___
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo