Geoapplications development http://rgeo.wikience.org Higher School - - PowerPoint PPT Presentation

▶
geoapplications development http rgeo wikience org
SMART_READER_LITE
LIVE PREVIEW

Geoapplications development http://rgeo.wikience.org Higher School - - PowerPoint PPT Presentation

Dec 31, 2022 183 likes •767 views

Geoapplications development http://rgeo.wikience.org Higher School of Economics, Moscow, www.cs.hse.ru 2 Why do we use maps? Pictures from .., (

slide-1
SLIDE 1

Geoapplications development http://rgeo.wikience.org

Higher School of Economics, Moscow, www.cs.hse.ru

slide-2
SLIDE 2

Why do we use maps?

2

Pictures from Лебедева О.А., Картографические проекции (методическое пособие), Новосибирский учебно-методический центр по ГИС и ДЗ, Новосибирск, 2000

slide-3
SLIDE 3

Map projection (leads to distortions)

3

  • https://en.wikipedia.org/wiki/Geoid
slide-4
SLIDE 4

Datum: aligning spheroid (ellipsoid)

4

Geoid/Earth surface Spheroid 1 Spheroid 2

slide-5
SLIDE 5

Putting all together

5

slide-6
SLIDE 6

Notes on CRS definition

6

select distinct coord_ref_sys_kind from epsg_coordinatereferencesystem

slide-7
SLIDE 7

Notes on CRS definition (2)

7

# R code # Convert (lat, lon) in meters to (lat, lon) in degrees library(raster) library(rgdal) # (lon, lat) in UTM projection – units are meters utm_c = cbind(470115, 6322515) utm_sp = SpatialPoints(utm_c, proj4string=CRS("+init=epsg:32637")) utm_sp # convert to (lon, lat) in WGS84 – units are degress global_sp <- spTransform(utm_sp, CRS("+init=epsg:4326")) global_sp

slide-8
SLIDE 8

Ellipsoids

8

slide-9
SLIDE 9

Datum

9

  • https://en.wikipedia.org/wiki/Geodetic_datum

https://en.wikipedia.org/wiki/World_Geodetic_System

slide-10
SLIDE 10

Coordinate systems (CS)

10

slide-11
SLIDE 11

Spherical geographic coordinate system

11

° ° ° °

Can you port that on Java from C++?

http://www.codeproject.com/Arti cles/15659/Longitude-Latitude- String-Parser-and-Formatter

slide-12
SLIDE 12

Cartesian coordinate system

12

slide-13
SLIDE 13

Cylindrical Transverse Cylindrical Oblique Cylindrical Secant Cylindrical Conical Secant Conical Planar Secant Planar

Projections

http%3A%2F%2Fwww.nps.gov%2Fgis%2Fgps%2F04_datums_coordinatesystems_65.ppt

slide-14
SLIDE 14

Projections

14

  • http://www.nps.gov/gis/gps/04_datums_coordinatesystems_65.ppt
slide-15
SLIDE 15

Why a map projection matters?

15

𝑒 = 𝑦1 − 𝑦2 2 + 𝑧1 − 𝑧2 2

slide-16
SLIDE 16

Equirectangular projection

16

Link 1

slide-17
SLIDE 17

Equirectangular projection

17

  • Link 1

The plate carrée (French, for flat square), is the special case where φ1 is zero.

slide-18
SLIDE 18

UTM: Universal Transverse Mercator

18

  • http://www.geo.hunter.cuny.e

du/~jochen/GTECH201/Lectur es/Lec6concepts/Map%20coor dinate%20systems/UTM%20a nd%20UPS.htm

slide-19
SLIDE 19

UTM: Universal Transverse Mercator

19

Apply a custom Transverse Mercator projection to each strip (zone)

Pic 1 Pic2 Pic3

slide-20
SLIDE 20

UTM: Universal Transverse Mercator

20

  • Custom means you have a separate cylinder for each strip (zone)
  • The cylinder is aligned to cross the spheroid along two lines (in red)
  • Each line is 180 km apart from the central meridian of a zone

Pic 1 Pic2 Pic3 Pic4

slide-21
SLIDE 21

UTM: zones

21

  • →
  • International Date Line)
slide-22
SLIDE 22

UTM: zones (another picture)

22

  • →
  • International Date Line)
  • PIC
slide-23
SLIDE 23

UTM: N & S

23

slide-24
SLIDE 24

UTM: False Easting and False Northing

24

Link 1

False easting False northing North zones 500,000 m none South zones 500,000 m 10,000,000 m

Link 2

False northing – no negative y-coordinates for all locations in southern hemisphere zones (10,000,000 m is enough) False easting – no negative x-coordinates west of a zone's central meridian. Northern zones do not use false northing: their y-coordinates are naturally positive. 𝒛 𝒚 𝒛 𝒚 = 𝟒𝟏𝟏 𝟏𝟏𝟏 𝒛 = 𝟓 𝟏𝟏𝟏 𝟏𝟏𝟏 𝒚 = −𝟑𝟏𝟏 𝟏𝟏𝟏 without false easting

slide-25
SLIDE 25

UTM: distortion

25

  • Link 1
slide-26
SLIDE 26

UTM: best practices

26

  • →
  • →

Link 1

slide-27
SLIDE 27

UTM: a limitation

27

  • °

°

  • Link 1
slide-28
SLIDE 28

Universal Polar Stereographic system

28

  • Link 1
slide-29
SLIDE 29

Military Grid Reference System (MGRS)

29

slide-30
SLIDE 30

Military Grid Reference System (MGRS)

30

  • Link 1
slide-31
SLIDE 31

NASA WW Example: MGRS Grid

31

MGRS

slide-32
SLIDE 32

Sinusoidal projection: very large scenes

32

slide-33
SLIDE 33

Sinusoidal projection

33 http://modis-atmos.gsfc.nasa.gov/MOD04_L2/grids.html https://en.wikipedia.org/wiki/Sinusoidal_projection

slide-34
SLIDE 34

EPSG database: www.epsg.org

34

Collection of CRS, CS, datum, etc. EPSG Dataset v8.7 contains: 5821 coordinate reference systems (CRS) from EPSG:200 to EPSG:69036405 126 coordinate systems (CS), CS != CRS 678 datum 14 variants of prime meridians 4232 CRS transformations (directly from → to) 1829 revisions ….. 21 tables

slide-35
SLIDE 35

CRS Text Formats

35

Try GeoTools, GDAL http://docs.geotools.org/stable/userguide/library/referencing/crs.html http://www.gdal.org/gdalsrsinfo.html

GDAL (gdalsrsinfo):

  • PROJ.4
  • OGC WKT
  • XML (GML based)
  • ESRI WKT format
  • Mapinfo style CoordSys format
  • simplified, etc.

http://spatialreference.org/

slide-36
SLIDE 36

raro@ubuntu-pelligrini:/mnt/hgfs/RS_DATA/Landsat8/LC81790212015146-SC20150806075046$ gdalinfo ./LC81790212015146LGN00_sr_band1.tif Driver: GTiff/GeoTIFF Files: ./LC81790212015146LGN00_sr_band1.tif Size is 8191, 8271 Coordinate System is: PROJCS["WGS 84 / UTM zone 37N", 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["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",39], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","32637"]] Origin = (224385.000000000000000,6322515.000000000000000) Pixel Size = (30.000000000000000,-30.000000000000000) Metadata: AREA_OR_POINT=Area Band_1=band 1 surface reflectance Image Structure Metadata: INTERLEAVE=BAND Corner Coordinates: Upper Left ( 224385.000, 6322515.000) ( 34d27'55.58"E, 56d57'49.93"N) Lower Left ( 224385.000, 6074385.000) ( 34d43' 3.37"E, 54d44'27.51"N) Upper Right ( 470115.000, 6322515.000) ( 38d30'26.82"E, 57d 2'42.39"N) Lower Right ( 470115.000, 6074385.000) ( 38d32' 5.80"E, 54d48'56.61"N) Center ( 347250.000, 6198450.000) ( 36d33'23.03"E, 55d54'25.94"N) Band 1 Block=8191x1 Type=Int16, ColorInterp=Gray Description = band 1 surface reflectance

GDAL output in WKT for Landsat 8 Moscow scene don’t worry due to font size – we are looking closer on it just in several slides

36

http://www.geoapi.org/3.0/javadoc/org/o pengis/referencing/doc-files/WKT.html

slide-37
SLIDE 37

PROJCS["WGS 84 / UTM zone 37N", 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["Transverse_Mercator"], PARAMETER["latitude_of_origin",0], PARAMETER["central_meridian",39], PARAMETER["scale_factor",0.9996], PARAMETER["false_easting",500000], PARAMETER["false_northing",0], UNIT["metre",1,AUTHORITY["EPSG","9001"]], AUTHORITY["EPSG","32637"]]

EPSG codes interpretation

37

One of the EPSG goals is to give short codes for frequent combinations to exchange metadata, the below WKT can be replaced by EPSG:32637 <geographic cs> <projected cs> UTM 37 North <prime meridian>

slide-38
SLIDE 38

EPSG codes interpretation (2)

38

slide-39
SLIDE 39

PROJ.4

39

The "no_defs" item ensures that no defaults are read from the defaults files. Sometimes they cause surprising problems.

http://lists.osgeo.org/pipermail/mapserver-users/2003-November/046863.html

gdalinfo -proj4 LC81790212015146LGN00_sr_band1.tif +proj=utm +zone=37 +datum=WGS84 +units=m +no_defs

slide-40
SLIDE 40

Notes on WGS84

40

As you can see from the EPSG database, WGS84 may mean (in different context) different entities:

  • Coordinate reference system
  • Datum
  • Ellipsoid
slide-41
SLIDE 41

Ellipse (1) – ellipsoid (2) – spheroid (3) – oblate spheroid, ellipsoid of revolution, reference ellipsoid, Earth ellipsoid (4)

41

  • Wikipedia: Ellipse

Wikipedia: Spheroid

  • Wikipedia: Ellipsoid
slide-42
SLIDE 42

Ellipse (1) – ellipsoid (2) – spheroid (3) – oblate spheroid, ellipsoid of revolution, reference ellipsoid, Earth ellipsoid (4)

42

(2) (3) (4)

SPHEROID["WGS 84", 6378137, 298.257223563, AUTHORITY["EPSG","7030"]]

  • )
  • https://en.wikipedia.org/wiki/Earth_ellipsoid

https://en.wikipedia.org/wiki/Reference_ellipsoid

slide-43
SLIDE 43

Reproject (raster=R, targetCRS=CRS, options)

43

→ →

slide-44
SLIDE 44

Reprojection: the process

44 http://demo.geo-solutions.it/share/OptimizingRasterReprojection.pdf

slide-45
SLIDE 45

Reproject: visualization in web maps

45

slide-46
SLIDE 46

Reproject use case: scenes from different projections

46

slide-47
SLIDE 47

Reproject: mosaic spanning several UTM zones

47

http://www.nrcan.gc.ca/earth-sciences/land-surface-vegetation/land-cover/north-american-landcover/9152 https://en.wikipedia.org/wiki/Canadian_Arctic_Archipelago

Lambert Conformal Conic (LCC) Ellipsoid Model: WGS84 (original Landsat projection is UTM)

slide-48
SLIDE 48

Reproject: why important

48

http%3A%2F%2Fwww.library.yale.edu%2FMapColl%2Ffiles%2Fdocs%2FUsing_Raster_Data.doc

slide-49
SLIDE 49

Vocabulary

49

Abbreviation Full CRS Coordinate Reference System SRS Spatial reference system (the same as CRS) CS Coordinate system (CS != CRS or SRS) Datum Spheroid axes and its relative position to the Earth WKT Well-Known Text EPSG European Petroleum Survey Group OGC Open Geospatial Consortium WGS84 World Geodetic System (1984) GIS Geographic Information System GML Geography Markup Language UTM Universal Transverse Mercator

slide-50
SLIDE 50

Main readings

50

Google for it

slide-51
SLIDE 51

Main readings

51

http://grassbook.org/

slide-52
SLIDE 52

Additional readings

52

slide-53
SLIDE 53

Other readings

53

QGIS tutorial http://www.qgistutorials.com/ru/docs/working_with_projections.html Other http://www.geo.hunter.cuny.edu/~jochen/GTECH201/Lectures/Lec6conc epts/Map%20coordinate%20systems/UTM%20and%20UPS.htm http://geology.isu.edu/wapi/geostac/Field_Exercise/topomaps/utm.htm https://www.maptools.com/tutorials/grid_zone_details https://en.wikipedia.org/wiki/List_of_map_projections

slide-54
SLIDE 54

Practical lesson 02

54

Resources: Data: EPSG, http://www.epsg.org/ (needs registration) Tools: To view EPSG: PostgreSQL, http://www.postgresql.org/ (open EPSG) + PgAdmin III

  • r MS Access or MySQL
slide-55
SLIDE 55

Practical lesson 02

55

GeoTools Maven quick start: GeoTIFF plugin: CRS info: EPSG plugin for GeoTools (be sure to have):

http://docs.geotools.org/latest/userguide/tutorial/quickstart/maven.html http://docs.geotools.org/stable/userguide/library/coverage/geotiff.html http://docs.geotools.org/stable/userguide/library/referencing/crs.html http://docs.geotools.org/latest/userguide/library/referencing/hsql.html Use Landsat 8 GeoTIFF for Moscow from previous lesson

slide-56
SLIDE 56

Practical lesson 02

56

  • 1. Use Java to print out CRS of GeoTIFF and specified one
  • 2. Use GDAL (gdalinfo, gdalwarp) to reproject rasters
  • 3. Use QGIS to visualize data in different projections

Data: Landsat 8 scenes Next training (#3): a short quiz on coordinates and projections

slide-57
SLIDE 57

Practical lesson 02: checklist

57

  • 1. Install OpenLayers Plugin
  • 2. Import vector/raster data in QGIS
  • 3. UTM Zone shapefile
  • 4. Project properties/CRS dialog box
  • 5. Pay attention to coordinates (false easting/northing)
  • 6. Identify features tool
  • 7. gdalinfo output
  • 8. try gdalwarp reprojection
slide-58
SLIDE 58