1

I download a certain geographic location (lon,lat) with buffer radius around for MODIS products, such as MCD43A2/MCD43A4 from Google Earth Engine using ee, such as:

 lon = 10.452175
 lat = 51.079206
 product = 'MODIS/006/MCD43A4'
 sensor = product.split('/')[-1]
 radius = 2000
 res = 500
 start,end = (datetime(2004, 1, 1, 0, 0), datetime(2004, 12, 31, 23, 59))


geometry = ee.Geometry.Point(lon, lat).buffer(radius)
col      = ee.ImageCollection(product)
features = col.filterDate(start, end).getInfo()['features']

http = urllib3.PoolManager(num_pools=10,
                            cert_reqs='CERT_REQUIRED',
                            ca_certs=certifi.where())


for feature in features:
    url = ee.Image(feature['id']).clip(geometry).getDownloadURL()
    #request.urlretrieve(url,OutDir.format(site_id,sensor)+feature['id'].split('/')[-1]+'.zip')
    r = http.request('GET', url)
    with open(OutDir.format(site_id,sensor)+feature['id'].split('/')[-1]+'.zip', "wb") as outfile:
        outfile.write(r.data)

than I get, using gdalinfo 2004_01_06.Nadir_Reflectance_Band1.tif :

Driver: GTiff/GeoTIFF
Files: 2004_01_06.Nadir_Reflectance_Band1.tif
Size is 10, 10
Coordinate System is:
PROJCS["MODIS Sinusoidal",
    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["Sinusoidal"],
    PARAMETER["longitude_of_center",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]]]
Origin = (727864.277675088495016,5682067.155503926798701)
Pixel Size = (463.312716528000010,-463.312716527000021)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  727864.278, 5682067.156) ( 10d25'45.53"E, 51d16'10.18"N)
Lower Left  (  727864.278, 5677434.028) ( 10d25'11.63"E, 51d13'40.26"N)
Upper Right (  732497.405, 5682067.156) ( 10d29'44.52"E, 51d16'10.18"N)
Lower Right (  732497.405, 5677434.028) ( 10d29'10.41"E, 51d13'40.26"N)
Center      (  730180.841, 5679750.592) ( 10d27'28.01"E, 51d14'55.22"N)
Band 1 Block=10x10 Type=Int16, ColorInterp=Gray

This does not refer to the input I made: 'lon = 10.452175, lat = 51.079206' at least for latitude.

How can I get around this, retrieving the original (projection) data subset, not re-projected, with the center (lon,lat) in the center of the subset? Or which would be the correct transformation to e done?

1

I have made some progress, changing GEE outputed geo.parameters. See the processing chain below:

Input coordinate: 10.453 51.07916667

gdalinfo on GEE extracted:

  Driver: GTiff/GeoTIFF
  Files: 2000_06_22.BRDF_Albedo_ValidObs_Band1.tif
  Size is 10, 10
  Coordinate System is:
  PROJCS["MODIS Sinusoidal",
      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["Sinusoidal"],
      PARAMETER["longitude_of_center",0],
      PARAMETER["false_easting",0],
      PARAMETER["false_northing",0],
      UNIT["metre",1,
          AUTHORITY["EPSG","9001"]]]
  Origin = (727864.277675088495016,5682067.155503926798701)
  Pixel Size = (463.312716528000010,-463.312716527000021)
  Metadata:
    AREA_OR_POINT=Area
  Image Structure Metadata:
    INTERLEAVE=BAND
  Corner Coordinates:
  Upper Left  (  727864.278, 5682067.156) ( 10d25'45.53"E, 51d16'10.18"N)
  Lower Left  (  727864.278, 5677434.028) ( 10d25'11.63"E, 51d13'40.26"N)
  Upper Right (  732497.405, 5682067.156) ( 10d29'44.52"E, 51d16'10.18"N)
  Lower Right (  732497.405, 5677434.028) ( 10d29'10.41"E, 51d13'40.26"N)
  Center      (  730180.841, 5679750.592) ( 10d27'28.01"E, 51d14'55.22"N)
  Band 1 Block=10x10 Type=UInt16, ColorInterp=Grayi

Change the flattening.geo.key in idl and save in tmp.tif

  tiff = read_tiff('2000_06_22.BRDF_Albedo_ValidObs_Band1.tif', geo=geo)
  % Loaded DLM: TIFF.
  IDL> geo.GEOGINVFLATTENINGGEOKEY = 0
  IDL> write_tiff,'tmp.tif', tiff, geotiff = geo

you can also edit the geoinformation of the *.tif with gdal_edit.py

gdal_edit.py -a_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' tmp.tif

gdalwarp with the original MODIS projection parameters

  gdalwarp -s_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +a=6371007.181 +b=6371007.181 +units=m +no_defs' -t_srs EPSG:4326 -r near tmp.tif out.tif
  Creating output file that is 15P x 8L.
  Processing input file tmp.tif
  0...10...20...30...40...50...60...70...80...90...100 - done.

checking the output:

gdalinfo out.tif
Driver: GTiff/GeoTIFF
Files: out.tif
Size is 15, 8
Coordinate System is:
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"]]
Origin = (10.414531289201173,51.099999995476345)
Pixel Size = (0.004987355776781,-0.004987355776781)
Metadata:
  AREA_OR_POINT=Area
  TIFFTAG_DATETIME=2019:07:11 10:28:07
  TIFFTAG_DOCUMENTNAME=tmp.tif
  TIFFTAG_IMAGEDESCRIPTION=IDL TIFF file
  TIFFTAG_RESOLUTIONUNIT=2 (pixels/inch)
  TIFFTAG_SOFTWARE=IDL 8.5.1, Exelis Visual Information Solutions, Inc., a subsidiary of Harris Corporation.
  TIFFTAG_XRESOLUTION=100
  TIFFTAG_YRESOLUTION=100
Image Structure Metadata:
  INTERLEAVE=BAND
Corner Coordinates:
Upper Left  (  10.4145313,  51.1000000) ( 10d24'52.31"E, 51d 6' 0.00"N)
Lower Left  (  10.4145313,  51.0601011) ( 10d24'52.31"E, 51d 3'36.36"N)
Upper Right (  10.4893416,  51.1000000) ( 10d29'21.63"E, 51d 6' 0.00"N)
Lower Right (  10.4893416,  51.0601011) ( 10d29'21.63"E, 51d 3'36.36"N)
Center      (  10.4519365,  51.0800506) ( 10d27' 6.97"E, 51d 4'48.18"N)
Band 1 Block=15x8 Type=Byte, ColorInterp=Gray

10.453 51.07916667

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.