1

R generated file:

enter image description here

Python generated file:

enter image description here

I want to georeference using external lat lng files in a MODIS data and then project it using R to WGS 84 coordinate system like how it happens using VRT file in Python. The code I could get so far shows wrong result in R when compared to result in Python.

    library(gdalUtils)   
    library(rgdal)
    fn<-'MOD11_L2.A2017080.0005.006.2017082010345.hdf'
    sds <- get_subdatasets('MOD11_L2.A2017080.0005.006.2017082010345.hdf')
    filename <- 'Path/file.tif'
    gdalinfo <- GDALinfo(fn, returnScaleOffset = FALSE)
    gdal_translate(sds[1], dst_dataset = filename)
    r <- raster(filename) #how to georefernce it using the lat long files in subdtaset of the hdf file
    wgs1984 <- CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0")
    projection(r) <- wgs1984
    extent(r)=c(xmn=130, xmx=170, ymn=-50, ymx=-20) # dont want to specify it since there are many files
    writeRaster(r, file="New_geo.tif", format='GTiff', overwrite=TRUE)
    plot(r)
    summary(r) 

#summary of file generated using R  
Min.    270.16  
1st Qu. 304.76  
Median  308.64  
3rd Qu. 310.92  
Max.    321.80  
NA's      0.00    
class       : RasterLayer   
dimensions  : 2030, 1354, 2748620  (nrow, ncol, ncell)  
resolution  : 0.0295421, 0.01477833  (x, y)  
extent      : 130, 170, -50, -20  (xmin, xmax, ymin, ymax)  
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0   
data source : Path\file.tif     
names       : file   
values      : 0, 65535  (min, max)  

#The files generated using python and warping shows this result. File generated in python gives this result   
        X1_LST  
Min.      0.00   
1st Qu.   0.00  
Median    0.00  
3rd Qu.   0.00  
Max.    321.58  
NA's      0.00     
class       : RasterLayer   
dimensions  : 1747, 2485, 4341295  (nrow, ncol, ncell)  
resolution  : 0.01217369, 0.01217369  (x, y)  
extent      : 133.8272, 164.0788, -43.75374, -22.48629  (xmin, xmax, ymin, ymax)  
coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0   
data source : Path\1_LST.tif   
names       : X1_LST   
  • Without either: a) your data or b) some data that shows the same effect or c) some summary info about your data its hard for us to help. You don't even explain how your result is "wrong". Please show as many outputs and summaries of your data objects as possible and indicate what seems "right" and what goes "wrong". – Spacedman Nov 24 '18 at 8:05
  • I have edited my question any more details required i will surely be able to tell. I think the problem here lies in the fact that I never warped the lat long file present in hdf file to the raster data set i only projected it to some projection. The python code does georeferencing using lat long vrt file. – ADIR Nov 24 '18 at 8:49
  • does python create the .vrt file, or do you get a vrt file? This can't be answered without seeing how that gets done. If you just get a vrt then raster() it. You can't apply GDAL's georeferencing logic for a raw coordinate data set in this way in R, but you could do it with gdalwarp directly – mdsumner Nov 24 '18 at 9:37
  • The python creates the vrt file. How can I do it in gdalwarping? lat=str(e)+'_lat.vrt' translate_lat='gdal_translate -of VRT %s %s' %(lat_str, latn) – ADIR Nov 24 '18 at 9:51
  • The code where the vrt is created in python is shown below. How can I use the lat long file in hdf to project it. – ADIR Nov 24 '18 at 10:03
0

Your R plot might fail because the data is swath data, i.e. the image is not parallel to meridians and latitudes. Using QGIS it displays like this:

enter image description here

The points are the GRINGPOINTS from the metadata.

So you can not simply take the min and max values for x and y, and squeeze the image pixels inside that frame.

You can use gdalwarp to batch-convert the MODIS data to georeferenced geotiffs as described in Batch conversion of Modis Swath (level 2) data and Batch for Mosaic of MODIS Swath

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.