texts:
I have the following script using rasterio:

if not os.path.isfile(out_dir + os.sep + name_adm1 + '.tif'):
    with rasterio.drivers():
        with rasterio.open(constants.ras_cmask) as src:
            b1 = src.read(1)

            profile = src.profile
            profile.update(
                dtype=rasterio.float64,
                count=1,
                compress='lzw')

        ndvi = b1.astype(float)
        ndvi[ndvi != int(num_ID)] = np.nan

        with rasterio.open(out_dir + os.sep + name_adm1 + '.tif', 'w', **profile) as dst:
            dst.write(ndvi.astype(rasterio.float64), 1)

However, when I run it, it produces a lot of messages like this:

03/13/2016 13:40:22  rasterio           DEBUG        Creating a not-responsible GDALEnv in drivers()
03/13/2016 13:40:22  rasterio           DEBUG        Got coordinate system
03/13/2016 13:40:22  rasterio           DEBUG        Got coordinate system
03/13/2016 13:40:22  rasterio           DEBUG        Input CRS: {'init': u'epsg:4326'}
03/13/2016 13:40:22  rasterio           DEBUG        Exported WKT: GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.0174532925199433,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]
03/13/2016 13:40:22  rasterio           DEBUG        Self CRS: {'init': u'epsg:4326'}
03/13/2016 13:40:22  rasterio           DEBUG        Got coordinate system
03/13/2016 13:40:22  rasterio           DEBUG        Got coordinate system
03/13/2016 13:40:23  rasterio           DEBUG        Option: ('BLOCKYSIZE', '1')

03/13/2016 13:40:23  rasterio           DEBUG        Option: ('COMPRESS', 'LZW')

03/13/2016 13:40:23  rasterio           DEBUG        Option: ('BLOCKXSIZE', '7200')

03/13/2016 13:40:23  rasterio           DEBUG        Option: ('TILED', 'FALSE')

Warning 6: driver ENVI does not support creation option BLOCKYSIZE
Warning 6: driver ENVI does not support creation option COMPRESS
Warning 6: driver ENVI does not support creation option BLOCKXSIZE
Warning 6: driver ENVI does not support creation option TILED
Warning 6: driver ENVI does not support creation option BLOCKYSIZE
Warning 6: driver ENVI does not support creation option COMPRESS
Warning 6: driver ENVI does not support creation option BLOCKXSIZE

Is there a way to either suppress these, or redirect to logging file?


texts: How to calculate azimuth from the south, latitude, departure with two points?

texts:
I'm making a Line class in python and the attributes we need to put in are the first point, end point, latitude, departure, distance, and azimuth (from the South). The points consists of eastings and northings. I understand that I need the azimuth first to calculate the latitude and departure. My question is how do I obtain the correct azimuth from the south with the two points?
I know the formula arctan((x2-x1)/(y2-y1)) must be used but that will not yield the correct azimuth from the south. How do I make the code without making another method?
Restriction: We can only use the math module.
I have this code so far.
class Line:
    def __init__(self,first,end):
        self.first = first #Type:Point (x1,y1)
        self.end = end #Type:Point (x2,y2)
        self.dist_x = self.end.x - self.first.x
        self.dist_y = self.end.y - self.first.y
        self.distance = math.sqrt((self.dist_x)**2 + (self.dist_y)**2)
        #self.azimuth = math.atan2((self.dist.x)/(self.dist.y))
        #self.latitude =
        #self.departure =


texts: How to edit a service definition file (SD) so as to change the Data source (User in oracle)

texts:
I am using Oracle 11g
ArcGIS 10.2

Considering I have an enterprise geodatabase and I am connected as user1/mypass@myIP/xe
I create a map  from my table "mylocations"
I  choose share as service-->Save service definition file --> No available connection and    choose the location that my sd is stored. (in arcmap 10.2)
Now my SD file knows where is my DB (ip) and from which instance (user1) and table (mylocations) my data comes from.
I have a development team where each user has his own schema in the same xe (user1Instance, user2Instance user5Instance etc)
Every user should be able to work with the Same sevice but with his OWN data.
As a result, I suppose I have to make 5 copies of the SD file and publish it in arcgis  for sever with different service names (service1, sevice2, ...sevice5)

The question can be raised:

How could I tweak in some way the different services (I suppose by editing in some way the SD files) so as each SD to know that it will retrieve its data from now on from a different location. For example SD1 should retrieve it data from table my location BUT from User1Instance, SD2 from user2Instance ..... sd5 form user5Instance and so on.

As a result when I publish these services, having registered the respective geodatabases, the server to be able to understand where each service should fetch its data.


texts: Expression to change an attribute each X number of rows

texts:
I'm using QGIS and I need an expression to change an attribute each X number of rows.
For instance I have the next attribute data:

I want to create an attribute (C) that assigns an ordered value each 3 rows. With the next result:


texts: Adding arrows to hydro line to show direction that river is flowing using ArcGIS Pro

texts:
I have a large number of rivers and streams on my map which do contain elevation data in the attribute table. I would like the lines on the map representing the rivers to have small reverse arrows sticking out of the line at an angle (i.e. like an arrow head) to show the direction that the rivers are flowing. I have tried to find a line symbology that could do this automatically for me but one does not appear to exist.
Is there any way around this problem?
I know that I could manually add hundreds of arrow heads on top of the hydro lines although I'm hoping that there is a better way.
I am using ArcGIS Pro with a Basic level license so I don't have all of the fancy tools.


texts: How to differentiate between land and sea in python?

texts:
I'm going to plot some positions of ships, and thus see which areas (out at sea) where there is no vessels, however, if I do this in an ordinary geographical area defined by a square, all the land area will be defined as an area without ships. This is ofc true, but how do I avoid it?

I have lots of lat/long coordinates of ship positions, and are primarily doing this in Python, but don't worry, a pseudocode will suffice.


texts: How can I read the TIFFTAG_DATETIME using GDAL in Python

texts: I am currently using GDAL from Python to read and transform coordinates and that works well. I would now like to be able to read some of the none geographic tags, for example, TIFFTAG_DATETIME using GDAL. I have searched the net but have been unable to find any examples. Can anyone help?


texts: How do I check Z-values in layer A with data specified in layer B in QGIS

texts:
I have information in a point layer with one Z-value per point. It is displayed in the screenshot as a label.
I have made an polygon layer with a "min Z value" per polygon. The polygons are displayed in the screenshot in brown, green and orange(ish). I have stored the "min Z value in the database of the polygonlayer.

I want to check if there are Z-values in the pointlayer that are less than specified in the "min Z value" column of the polygon layer.

Example: the "min Z value" for the orange(ish) polygon is 7.00;
I want to check if there are any values in the pointlayer within the perimeters of the orange polygon that are less then 7.00.


texts: Slope (Spatial Analyst) output shapefile having weird values

texts:
I'm trying to calculate the slope for a raster layer using ArcMap 10.2 and the Slope (Spatial Analyst) function. The raster I'm using is an NED 10m mosaic DEM, and the mosaic's datum is in WGS 1984. I am trying to calculate the degree slope, but my numbers are wonky. I know the output shapefile for the slope is incorrect because, according to these calculations, known flat areas have a slope of 71.9 - 89.9 degrees. My Z factor is 1 currently. I've tried reprojecting my DEM to NAD 1983 (2011) and change the Z factor, but I got the same results. Does anyone know what the problem I'm encountering is, and if so, how do I fix it?

One thing to note: I noticed my mosaic does not have an established Linear Unit (Layer Properties > Source > Spatial Reference (scroll)), which I found odd. Not sure if that ties into my problem or not, but I figured I'd mention it.


texts: Maximum number of returns per laser pulse in airborne LiDAR?

texts:
What is the maximum amount of returns a discrete aerial laser scanner is capable of recording?

For example in a LAS 1.2 metadata under Point Data Records there is the category
'Number of Returns' with minimum and maximum values listed:

  number_of_returns   1          6
  edge_of_flight_line 0          1
  scan_direction_flag 1          1
  classification      1         12
  scan_angle_rank   -30         29
  user_data           0          0
  point_source_ID    45         80
  gps_time 176457.257923 348033.538859
number of first returns:        1906575
number of intermediate returns: 2157
number of last returns:         1906575
number of single returns:       1881894
WARNING: there is 1 point with return number 6
overview over number of returns of given pulse: 1881894 45403 5441 593 76 6 0

As I got a warning I wanted to check what was the maximum number of returns the laser software and/or the processing software was capable of recording. Would the warning indicate that a return number would be erroneous?

I also looked at an example of a laser scanner (Riegl LMS-Q1560) and could not find any reference to maximum possible returns.


texts: Making black and clear raster image from polygon in QGIS?

texts:
I have made a polygon in QGIS (I just drew a polygon), but how do I make this into a raster image?

I want the raster image to be black (where the polygon was) and clear where the polygon wasn't.

I have tried using rasterize with a variety of options, but I always end up with an empty raster.


texts: ST_Intersection with only one polygon type

texts:
I am using this sql-query to execute an intersection:

SELECT ST_AsText(ST_Intersection (schablone.geom, eingangsdaten.geom)) --as geom2
     FROM schablone, eingangsdaten
     JOIN schablone b
     ON ST_Intersects (eingangsdaten.geom, b.geom)
     WHERE schablone.id = 4);

The result are different geometry types (polygon, multipolygon...). I need the result geometries with one geometry type. Is this possible?


texts: Why are JPG files invisible in Georeferencer?

texts:
I want to georeference a jpg file (map) with QGIS 1.8.0 on a Windows 7 (64 bit) machine.

The georeference window opens, but the map is not visible.
When converting the jpg file into a tif file (e.g. with gimp) everything works.
I tried it with several files. Always the same: The tif files are visible and can easily be georeferenced, while the jpg-files are not visible in the georeference window.

Any ideas?

Thanks for answers,
Ulli


texts: What will happen if we subtract two raster with different cell size using Raster Calculator?

texts:
I want to subtract a topography raster with a depth raster. The problem is the two rasters have different cell size. The topography raster have a square cell size of 5, 5 while the depth raster have a rectangular cell size of 10, 3.

How the Raster Calculator work with that way? Does it give accurate result?


texts: Difference between Feature Server and Map Server?

texts:
I have started to work with Esri maps.

What is the difference between Feature Server and Map Server?

(if possible with examples )


texts: Adding new map to current project using ArcGIS Pro .NET SDK?

texts:
The new ArcGIS Pro SDK for .NET is currently in beta for ArcGIS Pro 1.1, and I'm having trouble with map/project objects.

I can add layers to an existing map, but I can't seem to figure out how to add a new map to the project:

public async void AddMapToProject()
{
    //create new map and set basemap layer
    Map mymap = new Map("MapName");
    mymap.SetBasemapLayers(Basemap.Oceans);
    //there's also Map.Create("MapName"), not sure which to use

    //access current project
    Project myproject = Project.Current;
    //this doesn't work
    await myproject.AddAsync(mymap);

    //project has AddAsync(Item item) method, but I don't think
    //this is for adding maps as it doesn't accept a map object and
    //ArcGIS.Desktop.Core.Item.Create(map.uri) doesn't work
}

Any ideas? Can't seem to find an example or the correct method anywhere. The API reference for the SDK is here but doesn't offer many examples. Couldn't find much in the github samples / community samples either.


texts:
I have 5x5km grid rasters of monthly average wind speed for the UK. I want to convert these into yearly averages, what is the best way of doing so?

The only way I'm aware of is to polygonise the 12 rasters and export geometry/attribute data through mmqgis but that generates point data rather than polygon data. I've read in this which might do the trick: Converting point sets to polygon boundaries?

Is there a better/more efficient way around this?


texts: Import multiple KML files into single table with preserving layer name as (additional) column

texts:
I have bunch (1000+) of KML files like this:

<Document>
  <name>Shape 10</name>
  ...
  <Placemark>
    <name></name>
    <MultiGeometry>
     ...
   </MultiGeometry>
  </Placemark>
</Document>

which need to be imported into PostgreSQL database.

I can import all files using ogr2ogr:

find data -name '*.kml' | \
  xargs -n1 ogr2ogr -f "PostgreSQL" -dim 2 -nln shapes -append PG:dbname=${db}

but:

importing w/o -nln options creates one table per KML file with name from /Document/name element ("Shape 10") which is not ideal
importing w/ -nln creates single table but shape name ("Shape 10") is lost (not recorded in table - column name is empty probably because Document/Placemark/name is empty

How I could import 1000's of KML files into single table with preserving / storing layer name in output table?


texts:
I have an OVF layer defined as follows:

<OGRVRTDataSource>
        <OGRVRTLayer name="UserPoint">
         <SrcDataSource>ODBC:un/pw@CLO</SrcDataSource>
         <SrcSQL>select PKey as FID, Feature from PointData</SrcSQL>
         <GeometryField encoding="WKB" field="Feature"/>
  <FID>FID</FID>
        </OGRVRTLayer>
</OGRVRTDataSource>

The ODBC connection is to a SQL Express 2008 database. The PointData table is defined as follows:

CREATE TABLE [dbo].[PointData](
 [PKey] [int] IDENTITY(1,1) NOT NULL,
 [Feature] [geometry] NULL,
 CONSTRAINT [PK_PointData] PRIMARY KEY CLUSTERED
(
 [PKey] ASC
)WITH (PAD_INDEX  = OFF, STATISTICS_NORECOMPUTE  = OFF, IGNORE_DUP_KEY = OFF, ALLOW_ROW_LOCKS  = ON, ALLOW_PAGE_LOCKS  = ON) ON [PRIMARY]
) ON [PRIMARY]

This query

select *, Feature.STAsText() from pointdata

yields

1 0x00000000010C0000000000C062400000000000C06240 POINT (150 150)

"ogrinfo test.ovf --debug on" runs for a very long time returning many lines of this sort of stuff:

ERROR 4: Update access not supported for VRT datasources.
OGR_ODBC: EstablishSession(DSN:"CLO", userid:"un", password:"pw")
ODBC: SQLConnect(CLO)
ODBC: CatalogNameL: (null)
Schema name: (null)

OGR_ODBC: Table PointData has no identified FID column.
OGR_ODBC: Table UserMapData has no identified FID column.
OGR_ODBC: Table CHECK_CONSTRAINTS has no identified FID column.

before finally saying:

OGR: OGROpen(ODBC:un/pw@CLO/02C79D10) succeeded as ODBC.
ODBC: ExecuteSQL(select PKey as FID, Feature from PointData) called.
OGR_ODBC: Table SELECT has no identified FID column.
OGR: OGROpen(test.ovf/00E17090) succeeded as VRT.
Had to open data source read-only.
INFO: Open of `test.ovf'
      using driver `VRT' successful.
OGR: GetLayerCount() = 1

1: UserPoint
OGR: ReleaseDataSource(ODBC:un/pw@CLO/02C79D10) dereferenced and now destroying.
ODBC: SQLDisconnect()

GDAL version is 1.6.0.

I have tried using no FID tag and using SrcLayer instead of SrcSQL. Anybody know what I could do to make it see the FID? I'm guessing it would be a lot faster if it didn't check every table in the database for an FID...

Thanks in advance for your assistance.


texts: ProjNet Datum Transforms from EPSG:3785 to WGS84 y-error

texts:
I'm trying to write a method (in VB.NET) to transform datum from EPSG:3785 to WGS84. I'm currently making use of ProjNet in my solution so I used it to perform this transformation. I thought everything was going swimmingly at first; I plugged in my data and got my transformed data points out. When I checked the values the results were not exactly what I hoped.

To check the results I used QGIS to compare the data sets (before and after the ProjNet transformation) and noticed that there is a significant error term appearing in the y-axis.

.

In terms of what my code does: firstly there isn't a lot of documentation for Proj.Net but there was some helpful posts I used here:

How to perform transformations using ProjNet?
https://projnet.codeplex.com/wikipage?title=CreateProjection&referringTitle=FAQ

Similar to the other post, my call to the ProjNet transformation looks like so:

    Dim csf = new ProjNet.CoordinateSystems.CoordinateSystemFactory()
    Dim ctf = new ProjNet.CoordinateSystems.Transformations.CoordinateTransformationFactory()

    Dim sys3785 as ICoordinateSystem = csf.CreateFromWkt(_EPSG3785)
    Dim sys84 as ICoordinateSystem = csf.CreateFromWkt(_WGS84)

    Dim transformToWgs84 As ICoordinateTransformation = ctf.CreateFromCoordinateSystems(sys3785, sys84)

    Dim epsg3785Points = New Double() {_x, _y}
    Dim wgs84Points As Double() = transformToWgs84.MathTransform.Transform(epsg3785Points)

Attributes _x and _y are data points in EPSG3785. _EPSG3785 and _WGS84 are constants defined from the spatialreference.org resources as:

 Private Const _EPSG3785 = "PROJCS[""Popular Visualisation CRS / Mercator"", GEOGCS[""Popular Visualisation CRS"", DATUM[""D_Popular_Visualisation_Datum"", SPHEROID[""Popular_Visualisation_Sphere"",6378137,0]], PRIMEM[""Greenwich"",0], UNIT[""Degree"",0.017453292519943295]],PROJECTION[""Mercator""],PARAMETER[""central_me idian"",0], PARAMETER[""scale_factor"",1],PARAMETER[""false_easting"",0], PARAMETER[""false_northing"",0], PARAMETER[""latitude_of_origin"",0], UNIT[""Meter"",1]]"
 Private Const _WGS84 = "GEOGCS[""GCS_WGS_1984"",DATUM[""D_WGS_1984"",SPHEROID[""WGS_1984"",6378137,298.257223563]],PRIMEM[""Greenwich"",0],UNIT[""Degree"",0.0174532925199433]]"

Now ProjNet was throwing an error at first with the spacial reference WKT because EPSG:3785 didn't have defined "latitude_of_origin". I'm not too sure why so I appended it with a value of 0 as shown above. That's pretty much as far as I got with it before I got this transformation error.

Out of interest I decided to write a little test to see how the y-values differ between QGIS and the ProjNet transformations. I found that the difference across the y-axis was fairly linear the further the point moves from the origin (x-axis was always near exact).
.

Now I could write a little correction factor to the y-axis to re-align the data but that doesn't seem to solve the issue here. To me the problem exists either as an error in the projection strings or within ProjNet's transformation.

Is there anything apparent as to why this transformation error is appearing?


texts: Possible with different Time frame sizes in Time Manager?

texts:
I have 5 polygon shapefiles from 2019, 2020, 2040, 2060 and 2080 (predicting a certain development). t_start is DateTime datatype (2019-01-01, 2020-01-02 etc.)

It works perfectly if I use the 4 layers for 2020-2080 and set the Time frame size to 20 years.

However, when I add the 2019-layer and keep Time frame size 20 years, I go 2019-2039-2059-2079-2099 (fair enough). And if I choose 1 year Time frame size it shows every single years from 2019 to 2080 (not very much happening 2021, 2022,...2039, 2041, 2042,... 2059 etc.)

I have been playing around with end time and accumulate, but without any luck so far.

I want to show 2019, 2020, 2040, 2060 and 2080 - combine Time frame size 1 year and 20 years. Is it somehow possible to set "unequal" time frames?

I'm using QGIS 3.12.0 and Time Manager 3.4.


texts: Moving points to nearest raster cell with values in ArcMap

texts:
I have two types of data, one that is a point vector (occurrences), and one that is a raster (environment). The environmental data has both cells with numerical values or NA. Both of the types of data have the wgs84 coordinate system.

Due to the coarse scale of the environmental data, some of the occurrences fall just outside of the environmental cells with values and instead are within the cells with NA.

I want to be able to shift the geographic location of the occurrences from the cells with NA to the nearest cell with numerical values.

I want to shift them geographically, and not just sample the value of the nearest cell, because I am uploading the data into another program that builds a model based on the association between where the occurrences are and the environmental data they fall within.

How do I do this in ArcMap?


texts: Points interpolation on a polygon geometry

texts:
I would like to know how to interpolate point data on a polygon geometry.
I have two PostGIS tables:

the first one contains hourly rain values evaluated in the nodes of a regular grid (specified in a column with POINT type geometry)
the second one contains other values associated with the triangles of an irregular grid (specified in a column with POLYGON type geometry)

I would like to interpolate the hourly rain values (first table) on the polygons of the irregular grid (second table); then, I would like to save the results in the second table.
Is there a way to do this using some PostGIS functions?
I may use the IDW method, but it is only an hypothesis.


texts: Displaying PDF files associated to feature objects location.

texts:
Are there any other methods allowing user to display PDF files associated to points feature object other than using Actions function available to QGIS Layer Properties?

I have more than one thousand points feature object associated to former scanned map and summary sheet. I would like to implement a mtehod that allows the user when clicking, using the mouse to the point object, the associated map and summary sheet display. Both scanned maps and summary sheets are in PDFs format.

In the point feature attributes I have my full path locating my PDFs files (i.e maps and summary sheet) in two differrent columns.

How could I implemente this within QGIS?

Thank


texts: Why does ArcGIS "forget" my specified default GDB for a map document (mxd)?

texts:
I'm using ArcGIS 10.0 on Windows 7 64 bit Ultimate.

I've noticed that sometimes ArcMap's geoprocessing tools as well as scripts seem to forget my setting in the mxd (confirmed via arcpy.env.workspace and arcpy.env.scratchWorkspace.

Is there any way to get ArcMap to start using my specified GDB?

I'm familiar with the document properties and how I can change the default GDB for a map. The issue I'm seeing happens after I've set it and its been using an alternate GDB for a while.

I find that at random times, ArcGIS suddenly decides it wants to go back to using the "default default" (which annoyingly can never be changed. I can close and reopen the document, set the alternate GDB I use as the default againvia the properties dialog, but it will not move back to using the alternative GDB. Toolbox and scripts keep placing output in the default ArcGIS location (documents/ArcGIS/default.gdb).

I was however playing around with this and found manually settings the arc.env.workspace in the Python window seemed to magically bring it back (even though I was setting the same value?). Executing the verbatim sa tool I had just ran put output into the correct GDB even though the dialog settings nor the environment settings had "changed".

Has anyone else come across this behavior?


texts:
I have a feature
 Class with the columns KOMMUNE and GKR and some other columns.. Both columns are string. In the column KOMMUNE I have the value 0723 and 0724. In column GKR I have the value 0107 and 0108.
My problem is: If column KOMMUNE has the value 0723 i will add 2000 to the value 0107 in the colum GKR so the value becomes 2107 and if the value in KOMMUNE is 0724 I will add 3000 to the value 0108 in the colum GKR so the value becomes 3108.

Have tried with arcpy.CalculateField_management but I can not select the values with it.

import sys, string, os, arcpy

from arcpy import env

env.workspace = r"C:\Tom\ArcPyTest\Slå sammen kommuner"

#Kommune_1 = arcpy.GetParameterAsText(0)

#Kommune_2 = arcpy.GetParameterAsText(1)

#Nytt_Kommune_Nr = arcpy.GetParameterAsText(2)

Kommune_1 = 723

Kommune_2 = 720

Nytt_Kommune_Nr = "0722"

FcKommune = r"C:\Tom\ArcPyTest\Slå sammen kommuner\07_Vestfold inn.gdb\gkr2013_brutto07"

arcpy.env.overwriteOutput = True

arcpy.AddField_management(FcKommune, "KOM", "TEXT", '', '', 4)
arcpy.AddField_management(FcKommune, "GKR", "TEXT", '', '', 4)

arcpy.AddField_management(FcKommune, "KOMMUNE", "TEXT", '', '', 4)
arcpy.CalculateField_management(FcKommune, "KOMMUNE",  "!" + "KOMM" + "!", "PYTHON")

rows = arcpy.da.UpdateCursor(FcKommune, ["KOMMUNE", "GRUNNKRETS", "KOM", "GKR"])

for row in rows:
    if (row[0]) == str(Kommune_1):
        row[2] = str(Nytt_Kommune_Nr)
    if (row[0]) == str(Kommune_1):
        row[3] = (row[1] [4:8])
    if (row[0]) == str(Kommune_2):
        row[2] = str(Nytt_Kommune_Nr)
    if (row[0]) == str(Kommune_2):
        row[3] = (row[1] [4:8])

    if (row[0]) == str(Kommune_1):
        arcpy.CalculateField_management(FcKommune, "GKR", "[GKR] +2000" , "VB", "")
    if (row[0]) == str(Kommune_2):
        arcpy.CalculateField_management(FcKommune, "GKR", "[GKR] +3000" , "VB", "")

    rows.updateRow(row)


texts: Find the maximum value for a point among several overlapping vector layers in QGIS3

texts:
Dataset 1: around 80 vectorlayers with polygons. Each polygon has a value for a noise measurement. Some of the layers (up to ten at at time) are intersecting at some locations.

Dataset 2: a vectorlayer with sever hundred points.

The goal:
I want to add a column to the point dataset with the maximum value at this place from any of the polygons from the first dataset, as in: "look at the point in question, then look if any of the other layers have a polygon there. Then look at all the values from the noise column of all these polygons and write the highest one into a new column in the attribute set of the point layer".

The problem:
If I just merge the noise-layers I later can't use join_attributes_by_location with the point layer as I get a fix geometry error (I guess because the features are overlapping).
Do I have to split all the features which overlap with a feature in a different layer by all the features with which they overlap and then merge all layers and then merge all features with the same extent preserving only the maximum noise value for each shape, and then use that layer for join_attributes_by_location?
This sounds really complex and also as if it will take forever.

How can I achieve the situation described above (under goal)? Is there another spatial join I can use?

Note 1: I use QGIS 3 and know some Python and very basic R.
Note 2: I don't want to find the maximum value in a column in one attribute table, just to clarify.

Note 2: I'm new to SE, this is my first question. I already searched for similar problems but couldn't find anything. Forgive me if I overlooked something or forgot some piece of information.


texts:
I am trying to calculate curvature in QGIS using the tool of SAGA slope, aspect and curvature. I am working with DEM and I want to do some terrain analysis. However, the result of curvature layers are too bad. Below there is an image of the general curvature. Profile and plan give the same result.

Why this curvature is so spottiness and do you know what is going wrong?


texts: Access an Image from Image Collection in Google Earth Engine

texts:
I have created code which accesses  Sentinel 2 images, clip them according to a shapefile I have and calculate NDVI for the whole image collection.
Now, I want to visualize the images I have created, to see them as a layer and maybe after even export them.
 I tried to change the ImageCollection into a list and then use "get" in order to create variable out of the images, but it didn't work and I get this error message anytime I try to visualize the images:

 Cannot add an object of type <ComputedObject> to the map.

I tried it also on the image collection before all the calculations, when it was a bunch of Sentinel images, but I got the same error.

This is my code:

/**
 * Function to mask clouds using the Sentinel-2 QA band
 * @param {ee.Image} image Sentinel-2 image
 * @return {ee.Image} cloud masked Sentinel-2 image
 */
function maskS2clouds(image) {
  var qa = image.select('QA60');

  // Bits 10 and 11 are clouds and cirrus, respectively.
  var cloudBitMask = 1 << 10;
  var cirrusBitMask = 1 << 11;

  // Both flags should be set to zero, indicating clear conditions.
  var mask = qa.bitwiseAnd(cloudBitMask).eq(0)
      .and(qa.bitwiseAnd(cirrusBitMask).eq(0));

  return image.updateMask(mask).divide(10000)
  .copyProperties(image, ['system:time_start']);
}

// Map the function over one year of data
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2015-06-23', '2019-07-02')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .select('B2','B3','B4','B8','QA60')
                  .filterBounds(geometry)
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};
var clippedCol=dataset.map(function(im){
   return im.clip(geometry);
});

print(clippedCol);

//function to calculate NDVI
var addNDVI = function(image) {
  var ndvi = image.normalizedDifference(['B8', 'B4'])
  .rename('NDVI')
  .copyProperties(image,['system:time_start']);
  return image.addBands(ndvi);

};

//NDVI to the clipped image collection
var withNDVI = clippedCol.map(addNDVI).select('NDVI');

var colorizedVis = {
  min: 0.0,
  max: 1.0,
  palette: [
    'FFFFFF', 'CE7E45', 'DF923D', 'F1B555', 'FCD163', '99B718', '74A901',
    '66A000', '529400', '3E8601', '207401', '056201', '004C00', '023B01',
    '012E01', '011D01', '011301'
  ],
};

var listOfImages = withNDVI.toList(withNDVI.size());
print('List:',listOfImages);
var img1 = listOfImages.get(0);

Map.addLayer(img1,colorizedVis,'img1');

The end goal is to be able to visualize any image I want from the ImageCollection and to be able to export it.


texts: Python script to check if feature class exists in geodatabase if not convert feature class to feature class

texts:
I would like to check in one geodatabase and list the feature classes that are in it; check another geodatabsae for the list of feature classes, if they do not exist in the 2nd geodatabase then execute feature class to feature class.  If one does then just skip.  How can this be done with Python?  Here is what I have so far.  After I get my list how can I use arcpy.Exists to check another .gdb from my list, and if an element in the list does appear in that .gdb then skip it, however if it does not, then execute the function?

>>> import arcpy
... from arcpy import env
... import os
...
... # Set the workspace for the ListFeatureClass function
... #
... env.workspace = r'C:\GIS_Data\MXDs\ServerMxdsR\working\gis_dev01.gdb'
...
... # Use the ListFeatureClasses function to return a list of
... #  shapefiles.
... #
... fcList = arcpy.ListFeatureClasses()


texts: Link to another table in CartoDB infowindow

texts:
I have two layers in CartoDB, one with spatial data, school_count, and one without, heisman_winners.  I want to visualize my spatial layer but when a user clicks on a feature I want it to display data from my non-spatial table.  I have a common column between the two tables: 'school'.  On my non-spatial tables I have multiple records that correspond to each 'school' (i.e. 'class', 'name', 'position', 'percentage', 'points', 'year').  Can I set  a SQL query that runs when a user clicks a feature that will return each corresponding record from another table?  So when they click it will look like this:

'name' , 'position', 'points', 'percentage', 'year', 'class'
'name' , 'position', 'points', 'percentage', 'year', 'class'
'name' , 'position', 'points', 'percentage', 'year', 'class'
etc...

I fond this block but I am not sure exactly how to tweak it correctly.


texts: QGIS 3.16 offline work for search and rescue missions

texts:
I want to use QGIS 3.16 for search and rescue missions. I'm looking for a way to work offline. I've already tried XYZ, but every attempt fails after several hours. I need very large layers from the entire state of Rhineland-Palatinate, as we can be alerted throughout the entire state.
Can you give me a good tip or show me a way that I can work offline?


texts: Calculating overlap area that belong to each point in QGIS

texts:
I want to calculate the overlap area of buffered points in QGIS.
I have tried 'Polygon self-intersection' of the buffered points, but there is only 1 overlap area for each overlap.
For example in the image below, the selected area actually belongs to no 76 & 77, but there is only 1 polygon with ID 76|77.
So, my question is how to calculate the overlap area of each point?
For example point 76 have xxx m2 overlap area which is from overlap with points 77, 102, 75, and 51. Since the ID for each overlap area is unique, I find it is difficult to aggregate the overlap area that belongs to each point.


texts: Problem displaying OpenLayers KML layer in IE

texts:
I am serving the following files from a the same directory on localhost:

map.html
OpenLayers.js
output.kml

map.html simply creates an OpenLayers map and displays the KML file as per the OpenLayers example. I access it via http://localhost:8080/map.html and it all works fine in Firefox and Chrome but not in IE. In IE I get the navigation buttons but no map. The OpenLayers example does work in IE...

The contents of map.html are shown below:

<title>Map</title>

Any ideas?


texts: Delete numerous NoData rasters from database

texts: I've got >1000 rasters in a folder with about 10% containing no data. I'm looking for a method, preferably through ArcGIS (also open to QGIS, GDAL, etc.), to run through the folder and delete all of the no data rasters. Is there a script to take care of this? Logically, I could just delete all of the files under a certain size, but logic doesn't always apply to ArcGIS.


texts: How to extract POINT from a [LINE|MULTILINE]STRING

texts:
I'd like to create a vectorLayer in QGIS to display all the points that form a LINESTRING or a MULTILINESTRING stored in a PostGIS DB.

I think that I need to transform all points of LINESTRING (or MULTILINESTRING) to POINT.

First Question:

Is there any PostGIS function to do that?

Second Question:

If I want to exclude duplicate points from the result POINT table, how can I do that?


texts: Meaning of r.viewshed output numbers in QGIS?

texts:
I have completed a r.viewshed analysis of a point in GRASS7. I have imported the arising Geotiff into QGIS and when viewing the layer in the layer window (after expanding the +) I am presented with a number range= Min 88, Max 93.

Can someone please let me know what these numbers relate to?


texts: Classification validation: Trouble with the Confusion Matrix

texts:
I'm supposed to validate my supervised classification with a Confusion Matrix, but I'm having a hard time getting what I should use as a comparison.

At first, I used the same ROIs I created for classification to generate random samples throughout the image and then I reviewed it, making sure the pixels actually belonged to the different classes (that's what I called my "truth"). Thing is, my teacher said that wasn't ok and that there are known methods of obtaining your "truth". I tried to find it, but I just don't know the correct keywords to use on my searches. "Confusion Matrix" returns generic results explaining how to read it.


texts: Geoserver Error: No Such Feature Type in WFS-T

texts:
I was able to successfully create a Wfs_T app using the tutorial from https://medium.com/@goldrydigital/wfs-t-with-openlayers-3-16-6fb6a820ac58. I then decided to use this code to edit a file I'd made by creating a shapefile in QGIS, uploading into PostGIS, and then bringing it into GeoServer. To do this, all I did was change the code from

$.ajax('http://localhost:8080/geoserver/his_workspace/ows?', {
            type: 'GET',
            data: {
                service: 'WFS',
                version: '1.1.0',
                request: 'GetFeature',
                typename: 'his_layer_name',
                srsname: 'EPSG:3857',
                bbox: extent.join(',') + ',EPSG:3857'
            }
            })
var formatGML = new ol.format.GML({
            featureNS: 'his_workspace',
            featureType: 'his_layer_name',
            srsName: 'EPSG:3857'
            });

to

$.ajax('http://localhost:8080/geoserver/my_workspace/ows?', {
            type: 'GET',
            data: {
                service: 'WFS',
                version: '1.1.0',
                request: 'GetFeature',
                typename: 'my_layer_name',
                srsname: 'EPSG:3857',
                bbox: extent.join(',') + ',EPSG:3857'
            }
            })
var formatGML = new ol.format.GML({
            featureNS: 'my_workspace',
            featureType: 'my_layer_name',
            srsName: 'EPSG:3857'
            });

And I also changed

 $.ajax('http://localhost:8080/geoserver/his_workspace/ows', {
    type: 'POST',
    dataType: 'xml',
    processData: false,
    contentType: 'text/xml',
    data: payload
}).done(function() {
    sourceWFS.clear();
});

to

 $.ajax('http://localhost:8080/geoserver/my_workspace/ows', {
    type: 'POST',
    dataType: 'xml',
    processData: false,
    contentType: 'text/xml',
    data: payload
}).done(function() {
    sourceWFS.clear();
});

And Geoserver sends me the message:

org.geoserver.wfs.WFSException: No such feature type my_workspace:my_layer_name.

My describe feature type result is this:

<xsd:schema xmlns:my_workspace="my_workspace" xmlns:gml="http://www.opengis.net/gml" xmlns:xsd="http://www.w3.org/2001/XMLSchema" elementFormDefault="qualified" targetNamespace="www.my_workspace.com">
  <xsd:import namespace="http://www.opengis.net/gml" schemaLocation="http://localhost:8080/geoserver/schemas/gml/3.1.1/base/gml.xsd"/>
  <xsd:complexType name="my_layer_nameType">
    <xsd: complexContent>
      <xsd: extension base="gml:AbstractFeatureType">
        <xsd:sequence>
          <xsd:element maxOccurs="1" minOccurs="1" name="id" nillable="false" type="xsd:long"/>
          <xsd:element maxOccurs="1" minOccurs="1" name="geom" nillable="true" type="gml:PointPropertyType"/>
        </xsd:sequence>
      </xsd: extension>
    </xsd: complexContent>
  </xsd:complexType>
  <xsd:element name="my_layer_name" substitutionGroup="gml:_Feature" type="my_workspace:my_layer_nameType"/>
</xsd:schema>

My workspace, data store, and namespace all have the same name. The constraints and other settings that the tutorial mentions seem to be the same for my file. I tried changing the namespace to a URI and I still get the error message. I do wonder if it's because the geometry is in a point(geometry) type while in the tutorial it's simply in a geometry type.

Edit: here is the XML: Raw XML request: <Transaction xmlns="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd"><Insert><my_layer_name xmlns="my_workspace"><geometry><Point xmlns="http://www.opengis.net/gml" srsName="EPSG:3857"><pos srsDimension="2">-724011.5319171895 7284143.047464156</pos></Point></geometry></my_layer_name></Insert></Transaction>

Does anyone know why this isn't going through?


texts: osm2pgsql: database "gis" does not exist

texts:
I am currently attempting to use OSM Bright to quickly style some OSM data, so I can export tiles without spending hours writing Carto in TileMill.

This is my first time using osm2pgsql to import OSM (.pbf) data into PostgreSQL. It seems that I am running from one issue into the next (missing style file, no hstore, etc. ) and I have just hit a point where I am not sure how to progress anymore.

The following message keeps popping up:

Error: Connection to database failed: FATAL:  database "gis" does not exist

I have nowhere entered that I want to access a database called gis .

Executing the help command and looking at the docs I realize that:

 -d|--database The name of the PostgreSQL database to connect to (default: gis).

Which tells me that it is most likely looking for the default database, which I do not understand, as I have indicated to use my already-created database (ethiopia):

Here is the command I have entered:

F:\randomGISstuff\ethiopia>osm2pgsql -c -G -U postgres -S -hstore F:\randomGISstuff\ethiopia\default.style  -U postgres -d ethiopia ethiopia-latest.osm.pbf

If it is easier for you to look at a cmd screenshot, here it is:

The command that is featured on the OSM Bright page is definitely not sufficient in my case.

EDIT SOLUTION:

Because the solution is all hidden in comments, I will mention here what solved my issue:
1) Create a database and add the hstore to it:

CREATE EXTENSION hstore

2) make sure to set the password, not in the osm2pgsql statement, but before, like:

SET pgpassword=postgres

3) start osm2pgsql

4) Execute the following command:

 osm2pgsql -c -G -U postgres --hstore -S F:\randomGISstuff\ethiopia\default.style -d ethiopia F:\randomGISstuff\ethiopia\et hiopia-latest.osm.pbf


texts:
I am trying to scale a point layer in QGIS 3.10 around a certain point A. Ideally, after the transformation each point should be x metres closer to A, x being a fixed value and not a variable. All points in my layer are further than x metres from A.
It would be an affine transformation for each point individually, but not for the whole layer. I feel like I'm overlooking some very simple solution, but can't figure out how to go about this and couldn't find any answers on this site. I am open to R solutions as well.


texts:
I am trying to do a RBR (Relative Burn Ratio) and for this, I am using two Sentinel 2 L2A images (before and after the fire). The native CRS of the data is EPSG:32634 - WGS 84 / UTM zone 34N - Projected .

Since SNAP can work with archived Sentinel 2 imagery, I use the data as it was downloaded from Open Access Hub.

After I do the processing and the results are achieved (one RBR raster), I am exporting it in GeoTIFF format and trying to open in in QGIS. In QGIS, the layer does not line up with the basemap. So, I try to reproject the layer in SNAP and export it again... same problem.

For the workflow, I'm using QGIS 3.4 LTR but this does not matter that much, since I've also tried to display the results in ArcGIS 10.2.2 also. The basemap is Google Satellite with the CRS EPSG:3857 - WGS 84 / Pseudo-Mercator - Projected.

Reprojection didn't work(neither in ArcGIS, QGIS or SNAP) nor changing the CRS until the layer lined up with the basemap.

The discrepancy between the result and the basemap can be seen below.


texts: Creating new database in SpatiaLite and insert the data using Python in QGIS?

texts:
I want to create a new database in spatialite and want to connect to the database and execute the sql queries.

conn = db.connect('myDB.sqlite')
# creating a Cursor
cur = conn.cursor()
sql = 'select * from table_name'
cur.execute(sql)

Here i am connecting to the existing db. But i want to create a new db and connect to the db to execute the queries. How to do that using python.


texts:
I have lots of (multi)polygons generated with this tool: http://gis.ucla.edu/apps/click2shp/
It provides a zip with dbf, shp, shx, prj files. The prj file contains the following:

  GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137,298.257223563]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]

Which is EPSG:4326 as I understand.

Polygons are loaded to the Postgres + PostGIS database using some third party library (unknown to me). And it is accessible as WKB or WKT using SQL.
The polygons are spread all over the world - USA, Australia, Europe, potentially anywhere.

What I need is the projected coordinate system for each of these polygons to calculate their areas.

So far I've tried exploring
http://www.epsg-registry.org
that allows me to identify the CRS using some point or bounding box.
But maybe there's a better (preferably automated) solution.


texts: Does the tag bicycle=yes mean that its a preffered path for bicycles?

texts:
If a path is tagged with

bicycle=yes

Does it mean that its a path where bicycles are allowed or is it actually a path that is preferred for bicycles (due to some reason like bicycle infrastructure)?

Edit:
I understand that the tag implies access but wanted to know if the tag can be used as an indication of a road more preferable for biking compared to roads thats don't have this tag.

Edit: To clarify my question further

I was mapping bicycle ways for a city and was wondering, since not every path has a bicycle tag, the ones that do have bicycle=yes, would they be more suitable for biking compared to the ones that don't have any' bicycle=*` tag (except the ones that imply cycling isn't allowed).


texts: SpatialJoin with "closest" option ok in python script but not when in arcmap toolbox

texts:
I have the following problem. I want to create a toolbox to obtain the distance between a polygon (Bldg.shp) and a point (pts.shp). I have started with a simple following python script:

import arcpy

bldg=r"C:\0_prg_dvp\SpatialJoin\Bldg.shp"
Pts=r"C:\0_prg_dvp\SpatialJoin\pts.shp"
Output=r"C:\0_prg_dvp\SpatialJoin\Test"
arcpy.SpatialJoin_analysis(bldg, Pts, Output+".shp", "JOIN_ONE_TO_ONE", "#", "#", "CLOSEST", "#", "HorzD")

All works fine. If I try to modify it to work in a toolbox as follow:

 import arcpy

 #bldg = arcpy.GetParameterAsText(0)
 #Pts = arcpy.GetParameterAsText(1)
 #OutputShp = arcpy.GetParameterAsText(2)
 #Output="C:/0_prg_dvp/SpatialJoin/"+OutputShp

 bldg=r"C:\0_prg_dvp\SpatialJoin\Bldg.shp"
 Pts=r"C:\0_prg_dvp\SpatialJoin\pts.shp"
 Output=r"C:\0_prg_dvp\SpatialJoin\Test3"

 arcpy.SpatialJoin_analysis(bldg, Pts, Output+".shp", "JOIN_ONE_TO_ONE", "#", "#", "CLOSEST", "#", "HorzD")
 #arcpy.SpatialJoin_analysis(bldg, Pts, Output+".shp", "JOIN_ONE_TO_ONE", "#", "#")

Arcmap crashes without any error messages, the output file is created but is empty. If I remove the "Closest" option, the toolbox script works fine (but I don't have the distance I want)...

Does anyone experienced this and have a fix? Am I missing something?

Thanks for your help


texts: Remove time stamp from date field using field calculator

texts: I'm trying to copy the dates from a field called InspectionDate (date data type) to new field called CompletedDate (date data type) but I don't want the time stamp to be copied over.  How can I copy just the date in a mm,dd,yyyy format using field calculator.  For example copying 1/1/2016 12:00:00 AM to new field as 1/1/2016.


texts: QGIS: Using Split Vector Layer how can I name the new files based on the unique value?

texts: I'm using Natural Earth 50m data and I want separate country shp layers so I can export to Illustrator. I can classify by Long_Name easily. I can Split Vector Layer on Long_Name easily. What I can't do is name the resulting files based on their unique Long_Name value. What am I missing?


texts: How to change text legend orientation from right to left reading with arcpy

texts:
I working on mxd with ArcMap 10.3 . The default's legend properties is- Layout - left to right reading (Marked with a red ellipse) :

when it unmarked i get this legend:

i want to know if there a way to mark the default box automatic from right to left reading option so the result will be like that (Very important is that all the legend items will be in straight line from left side):

Does anyone know of a way to do this with arcpy?

UPDATE:

I tried MrBubbles answer and the result was:

while i wanted result like that (i made it manually for the example):

The answer of MrBubbles didn't help me because it very important that all the words in the legend will be right to left configuration
How can i do it?


texts: Why is QGIS holding references to layers that are not in layer list

texts:
I have a QGIS project, which has about 30 layers in the layer list. I removed one of the layers from the list, saved the project, closed QGIS.
(The layer itself was from a Geopackage database, which has several other layers in it, which are still in the layer list).
I have reopened the project (few days later) and - the layer is still not in the Layers list (as you would expect).
However - navigating to the Vector Geoprocessing tools - The layer is listed there as a potential input for the Geoprocessing tool (if i select it as an input, the Geoprocessing tool simply fails after about 20 seconds).
I can add the layer back in, and it then appears twice in the layer list of the geoprocessing tools.
I can't see why this layer would appear as an input layer to the Geoprocessing tools, when it is clearly not in the Layers list.
Im hoping for some help in identifying why. More specifically:

A python query on either a standard Geoprocessing tool (like Buffer) and/or the layer list to see if it is listed there
Console logging that I could implement on QGIS Project startup which would also help identify where it is being loaded as part of the project.
Any other ideas?


texts: Object Inspector Button not enabled in python console of qgis

texts: I am using QGIS2.6 version. I loaded an existing py file in a editor. Though I added a codes in the editor, the object inspector button was not enabled. Please help me in figuring out what would be the problem.


texts:
I'm manipulating an OSM dataset of Chicago to find the optimal path for heavy duty vehicles. Are there some roads restrictions? Which class is not allowed for heavy vehicles?

These are road classes of the OSM map:

"secondary"
"motorway"
"trunk_link"
"tertiary"
"unclassified"
"secondary_link"
"tertiary_link"
"trunk"
"primary"
"residential"
"primary_link"
"motorway_link"

I asked to get a dataset of road restriction here and I didn't get a pertinent answer.


texts:
i`m working with QGIS 3.10 and 3.12, trying to use the WhiteboxTools plugin. I followed the WBT-Manual (point 3.4 QGIS Plugin).The Repository details dialog box is mentioned as "https://plugins.bruy.me/plugins/plugins.xml" but it creates the Error:
"Error transferring https://plugins.bruy.me/plugins/plugins.xml?qgis=3.12 - server replied: Not Found"
Is there a new link needed?
My friend got the same problem. Worked 1 month ago, but now he has the same error.
I'm using Windows 10.


texts: Moving map content/ map signatures when zoom

texts:
Got a problem with my leaflet map especially when using it on mobile devices.

When using the zoom (zoom in and out) my signatures are moving out of the screen. after finishing the zoom the signatures get their correct location.

is this the normal behavior of leaflet cuz of the layers? and how to fix this?

when adding a L.Marker f.e. for GPS location or for search results the marker doesnt move but is fixed on his position while zoom. would like to have the same behavior for my signatures.

signature = point signatures --> points --> in this case POIs

version of leaflet : 0.7.5


texts: QGIS: Extract rows of grape vines from inter-row vegetation using a 3D point cloud

texts:
I have some multi spectral imagery (raster data) from a UAV of vineyards in Victoria, as well as corresponding 3D point cloud (also raster) of the vineyard topography, based off ultrasonic sensing, which picks up the height difference between the vines and the ground between the rows.

The article linked below discusses a method for extracting the vine pixels using the height differentials expressed in the point cloud.

http://www.mdpi.com/2072-4292/9/2/111/

Specifically they divide the resulting imagery up into a grid consisting of 10m x 10m squares and calculate a reference height for each square (using to lowest 10% of pixels in the point data). Pixels exceeding the reference height reflect the vine pixels and are extracted.

I'm looking for some guidance as to how I might do this using QGIS. I'd like to generate a binary mask that could be used to say, extract NDVI data for the vines only.

I'm fairly well versed with QIS and how to drive it, but this has me stumped.


texts: Using QGIS API to check which layers are being used in layout

texts:
I have a QGIS project with multiple layouts (in total about 30) in my layout manager. I used of the "Follow map theme" then "Lock layers" options in the Map item properties to make sure that each layout is displaying the correct layer. However after I am done making all the layouts I would like to make sure that they are displaying the correct layers.

Is there a way in QGIS to verify which layer is being displayed in the layout? Making all the layouts manually is prone to error. This option is absolutely necessary to make big map making projects possible.

I am also open to using the Python API to do this.

I have made some progress:

project = QgsProject.instance()
manager = project.layoutManager()
layouts = manager.layouts
layout = layouts[0]
items = layout.items()

#find map object
for item in items:
    if isinstance(item, qgis._core.QgsLayoutItemMap):
        map_object = item

#Here need to find map object somehow
layers = map_object.layers

The above will extract the layers as QgsVectorLayer or QgsRasterLayer objects. Unfortunately this is not enough, I would also like to extract the position of the QgsVectorLayer in the layer tree. I found out how to go down a layer tree from this post. However when I get down to the object of interest it returns it as a QgsLayerTreeLayer.

It is possible to extract the underlying layer from the tree layer by doing this:

layer_from_layer_tree = tree_layer.layer()

And it is apparently possible to also do this:

layer_from_layer_tree==layer_from_map_layout
>>> True

But I am not sure if this is the recommended way to compare two layer objects, and if objects which are not equal will return False.

In any case, I think it is very convoluted. Verifying which layers are in your layout should be a basic functionality available directly in the GUI.


texts: Get bounds of featureGroup with multiple features

texts:
I tried looking into a solution for this myself, or an alternative way to load my layers, but I'm simply unable to.

I first load a geoJSON file through L.GeoJSON.AJAX(), then add the features (polylines and markers) that I need to a L.featureGroup(), but I'm having issues using the featureGroup.getBounds() function when doing this.

How can I get the bounds of all the features added to a featureGroup like this?

Edit:
Added code. Completely slipped my mind when asking the question to include this.

function addMapData(id, map) {
  var newIcon = L.icon({
    iconUrl: './img/icons/ring.png',
    iconSize: [84, 80],
    iconAnchor: [5, 70],
    popupAnchor: [10, 10],
  });

  var mapData = L.featureGroup();

  new L.GeoJSON.AJAX('./data/fjelltrimmen/fjelltrimmen.vefsn.geojson', {
    local: true,
    style: {
      weight: 5,
      color: '#c0392b',
    },
    onEachFeature: function(feature, layer) {
      if (feature.properties.id === id) {
        // Polyline
        layer.addTo(mapData);

        // Map marker
        L.marker([
          feature.geometry.coordinates[feature.geometry.coordinates.length - 1][1],
          feature.geometry.coordinates[feature.geometry.coordinates.length - 1][0],
        ], { icon: newIcon }).addTo(mapData);
      }
    }
  });

  mapData.addTo(map);
  map.fitBounds(mapData.getBounds());
}


texts: TypeError: "unsupported operand type(s) for /: 'int' and 'list'" in QGIS?

texts:
In a academic project I've produced a small python script for watershed analysis. This script relates shapefiles attributes and generates a .txt file.

I am trying to adjust the script to use within the QGIS environment. When I use it, it shows the following error:

  unsupported operand type (s) for /: 'int' and 'list' See log for more details

In logs we can see:

2019-01-03T21:10:57 2   Uncaught error while executing algorithm
        Traceback (most recent call last):
          File "C:/PROGRA~1/QGIS2~1.18/apps/qgis-ltr/./python/plugins\processing\core\GeoAlgorithm.py", line 203, in execute
            self.processAlgorithm(progress)
          File "C:/PROGRA~1/QGIS2~1.18/apps/qgis-ltr/./python/plugins\processing\script\ScriptAlgorithm.py", line 381, in processAlgorithm
            exec((script), ns)
          File "<string>", line 123, in <module>
          File "<string>", line 35, in densidade_rios
        TypeError: unsupported operand type(s) for /: 'int' and 'list'

Here is the code:

from qgis.core import *
import math

# Function
def densidade_rios(valor_total_rios, valor_area_bacia):
    dr = valor_total_rios / valor_area_bacia
    return dr

# Loading dreinage and watershed layers
dlayer = QgsVectorLayer(Drenagem, "drenagem", "ogr")
blayer = QgsVectorLayer(Bacia, "bacia", "ogr")

# Dreinage attributes
dfeatures = dlayer.getFeatures()
# Listing values per column
listdfeatures = list(zip(*dfeatures))

# Max and min order values
maxordem = max(listdfeatures[0])
minordem = min(listdfeatures[0])

# Counting features of channel hierarchy (dreinage layer values)
ordem = minordem
hierarquia = []
while (ordem <= maxordem):
    contagem = listdfeatures[0].count(ordem)
    hierarquia.append(contagem)
    ordem += 1

# Total of rivers - quantity of first order rivers (dreinage layer values)
totalrios = hierarquia[0]

# Watershed attributes
areabh = blayer.getValues(blayer.fields()[0].name())[0]

# Using function
dr = densidade_rios(totalrios, areabh)

Attribute table values of drenagem layer:

Attribute table values of bacia layer:


texts: Downloading data from Google maps-based apps

texts: Is there any way to download customer’s data in Google maps based app? Like the polygons in this map - http://tromso-parkering.no/parkering/kartoversikt/.


texts: Calculating a bounding box by its center point and ratio, the map zoom and desired pixels width in Leaflet?

texts:
I'm trying to add dynamically an overlay on a leaflet map.

I know the

overlay center - {lat:50.8466,lng:4.3528}
overlay ratio - 16/9
map zoom - 15
desired pixels width - 200px

How can I calculate a bounding box based on those parameters?

EDIT

I ended up working with the map size (in pixels) rather that with the map zoom; with the method containerPointToLatLng():

  Given a pixel coordinate relative to the map container, returns the
  corresponding geographical coordinate (for the current zoom level).


texts: Who should be responsible for editing GIS data?

texts:
I work for a municipality with an enterprise GIS. What would be the argument for performing all data editing in the GIS department? Conversely, what would be the argument for data updates being handled at the point of entry (e.g. DPW: sewer & water; Assessing: parcels).

I am interested in the ideal case, regardless of staffing and money.


texts: Time Manager to animate time series in multiple columns?

texts:
I have a vector file of Africa containing all its countries. I have a CSV file containing a time series of monthly precipitation per African country. Is there a way to use Time Manager (or anything else) in computing a series of images containing all countries with their corresponding precipitation value over time?

So what I would like to achieve is an image of Africa with all countries having a different color based on a precipitation-value, and this for every month over a period of about 30 years. Obviously I would like to avoid having to make 12 x 30 separate images myself.

Tutorials that I find only show how to use Time Manager with a CSV that contains spatial data, which is not the case for me. If I add all the precip-data as columns to the attribute table of the Africa shapefile, I have the problem of only being able to select one column for Time Manager to run. If I do it the other way around I can only do it for one country in stead of all countries simultaneously.

Any suggestions as to how I can achieve this?


texts: How to insert attributes into a WFS using OpenLayers?

texts:
So I have a web app that's running and is able to add points and lines to existing Postgis tables by getting data from Geoserver. I need the following one specific bit of functionality added:

I want to allow users to also be able to add attribute information to each point on the map interface. As in, each point they draw on the map, there needs to be something that allows the user to also enter some text data to the column.

I've tried reading a few questions in here but none of them provide a solution for a point vector layer.

How to do this?

The bit I have for loading and posting the WFS point features is:

var vectorSource2 = new ol.source.Vector({
loader: function(extent, resolution, projection) {
var url2 =    'http://localhost:8080/geoserver/wfs?service=WFS&' +
'version=1.1.0&request=GetFeature&typename=BFTchambers:chamber2&' +
'outputFormat=text/javascript&format_options=callback:loadFeatures2' +
'&bbox=' + extent.join(',');
$.ajax({url: url2, dataType: 'jsonp', jsonp: false})
.done(function(response) {
    pointWFS = new ol.format.WFS(),
    sourceVector2.addFeatures(pointWFS.readFeatures(response))
    });
},
strategy: ol.loadingstrategy.tile(ol.tilegrid.createXYZ({
maxZoom: 20
})),
});

window.loadFeatures2 = function(response) {
console.log("Point features were drawn");
vectorSource2.addFeatures(geojsonFormat.readFeatures(response));
};
var formatWFS2 = new ol.format.WFS();
var pointGML = new ol.format.GML({
featureNS: 'http://geoserver.org/bftchamber',
featureType: 'chamber2',
});

var pointWFS = function(p,f) {
switch(p) {
case 'insert':
    node = formatWFS2.writeTransaction([f],null,null,pointGML);
    break;
case 'update':
    node = formatWFS2.writeTransaction(null,[f],null,pointGML);
    break;
case 'delete':
    node = formatWFS2.writeTransaction(null,null,[f],pointGML);
    break;
}
s = new XMLSerializer();
str = s.serializeToString(node);
$.ajax('http://localhost:8080/geoserver/wfs',{
    type: 'POST',
    dataType: 'xml',
    processData: false,
    contentType: 'text/xml',
    data: str
    }).done();
    console.log(" point features were posted to server");
}
case 'btnDrawPoint':
    interaction = new ol.interaction.Draw({
        type: 'Point',
        source: layerVector.getSource()
    });
    map.addInteraction(interaction);
    interaction.on('drawend', function(e) {
        pointWFS('insert',e.feature);
    });
    break;


texts: QGIS: ideas on how to show unsupported geometry data type

texts:
I am working on a university project and I need to visualize data from PostgreSQL. The data type is new, geometry with a time attached to it, which it is available because of extension in PostgreSQL.

I need to find a way to visualize this data in QGIS (create a plugin), and I couldn't find relevant information online. This blog post is similar to what I am trying to achieve, but using the data type mentioned before.
* https://anitagraser.com/2017/08/14/movement-data-in-gis-7-animated-trajectories-with-timemanager/

Does anyone here have ideas on how to approach?


texts: Changing style of different zones in OpenLayers?


texts: Bare World Satellite Map Base Layer and Tile Server?

texts:
I'm fairly new to mapping, and am working on a custom data visualization that will be presented on top of a Leaflet map. We'll be placing our own geoJSON polygons on the map (including labels) and would like to have a bare satellite map base layer as the background (i.e. no administrative or other labels).
I've found the following satellite map layers and tile servers:
Mapbox: https://www.mapbox.com/maps/satellite#the-data - requires registration, and offers up to 50,000 map / tile server requests per month in their free tier.
Maptiler: https://cloud.maptiler.com/maps/hybrid/ - requires registration and I believe offers up to 100,000 tile server requests per month. Fortuantely too they offer the 'bare' tiles here - https://cloud.maptiler.com/tiles/satellite/ - which is exactly what we want (they also offer custom maps  - which ironically, if using the above satellite tiles with all administrative labels removed via their advanced editor, is no longer in their free tier).
Esri ArcGIS World Imagery: https://www.arcgis.com/home/webmap/viewer.html?layers=10df2279f9684e4a9f6a7f08febac2a9 - which I believe is free to use via a developer account (not sure what the monthly restrictions are).
We really like the map and water color of the Maptiler and Mapbox maps.
Since we're looking for a bare base layer map with no administrative or other labels, is there any other option for including a satellite base layer like the ones from Maptiler or Mapbox in our Leaflet map that doesn't require registration with one of the above services?


texts: How to return output a raster layer in a custom processing script in Pyqgis

texts:
I've written a couple of custom scripts for vector shapefiles by following the template provided in QGIS. I am now trying to write a script which will use gdal algorithms to process layers. I have found a great guide here on how to adapt the template for raster layers. My question is if I use a gdal algorithm, how do I pass it to the output to be returned by the script?

For example, if I have rasterized = processing.run("gdal:rasterize", ... in the calculations part, how do I then pass that to the output_path_raster_a = self.parameterAsOutputLayer...


texts: Finding line segment between line intersections

texts:
Given 1 non-self-intersecting line (red) and 2 other non-self-intersecting lines (light-green), what is the best way to find the coordinates of the line in between the 2 intersections?
I threw this image together to better illustrate what I'm saying.
Preferably to be implemented in PHP code, does anyone have an algorithm to offer as to how to find the coordinates of the red line in between the 2 red/light-green line intersections?
Coordinates will either be in lat/lon or X/Y plane format though this shouldn't affect the algorithm.


texts: Count number of points within certain ranges from set of points using PostGIS?

texts:
I have the following, which gives me the number of customers within 10000 meters of any store location:

SELECT COUNT(*) as customer_count FROM customer_table c
WHERE EXISTS(
SELECT 1 FROM locations_table s
WHERE ST_Distance_Sphere(s.the_geom, c.the_geom) < 10000
)

What I need is for this query to return not only the number of customers within 10000 meters, but also the following. The number of customers within...

10,000 meters
more than 10,000, but less than 50,000
more than 50,000, but less than 10,0000
more than 100,000

...of any location.

I'm open to this working a couple of ways. For a given customer, only count them one time (the shortest distance to any store), which would count everyone exactly once. I realize this is probably pretty complex. I'm also open to having people be counted multiple times, which is really the accurate values anyway and think should be much simpler.


texts: problem combining geometries from all features in a layer as one geometry

texts:
I am trying to combine all features of a Polygon layer in QGIS in a python plugin, and am running into this error when I try to run the combine method:

"AttributeError: 'NoneType' object has no attribute 'combine'
Python version:
2.7.5 (default, May 15 2013, 22:44:16) [MSC v.1500 64 bit (AMD64)]

QGIS version:
2.10.1-Pisa Pisa, d20c5b7"

I've looked at this question before but am not sure what I am doing incorrectly:

from qgis.core import *
from PyQt4.QtCore import QSettings, QTranslator, qVersion, QCoreApplication
from PyQt4.QtGui import QAction, QIcon, QFileDialog
# Initialize Qt resources from file resources.py
import resources, processing
# Import the code for the dialog
from save_attributes_dialog import SaveAttributesDialog
import os.path

Here is my code for the run method:

def run(self):
    layers = self.iface.legendInterface().layers()
    layer_list = []
    for layer in layers:
        layer_list.append(layer.name())
        self.dlg.comboBox.addItems(layer_list)
    # show the dialog
    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed
    if result:

        selectedLayerIndex = self.dlg.comboBox.currentIndex()
        selectedLayer = layers[selectedLayerIndex]

        features = selectedLayer.getFeatures()
        geoms = QgsGeometry()

        for feature in features:
            geoms = geoms.combine(feature.geometry())

I ran the loop just printing out the areas of each feature before which worked fine, so I'm guessing it is able to access the features without a problem, but am not sure what to do about the NoneType issue.


texts: Correct order of nodes in Overpass query

texts:
I am struggling to get the correct order of nodes out of an Overpass query.

I am aware of this similar question How do you get the nodes of an area in Overpass API in the "right" order? but I am not sure if I understand the solution hinted there. I wish the text was a bit more explicit how this is solved in the end.

Let's for example use way ID 35565161:
http://www.openstreetmap.org/way/35565161#map=15/49.3949/4.0426
http://api.openstreetmap.org/api/0.6/way/35565161

What I found is that running my query "api/interpreter?data=way(35565161);>;out body;" does not return the nodes in the correct order.

I think I understand the answer to the question linked above that I have to first run a query on the way to get the order of nodes, then run another query on the node details to get the node coordinates.

With those two available iterate over the way, referencing the detail results from the second.

That seems a bit too much work - I guess there is an easier way to do this in one query?


texts: Best flow routing algorithm for determining location where spill enters stream?

texts:
What would be the best flow routing algorithm for determining the quickest (and by this I suppose I mean shortest) path that a spill would take when traveling to a stream?

My best educated guess would be that the D8 would provide the quickest path for a spill traveling from land to the stream.  It minimizes dispersion (which would most likely slow the spill down) and maximizes velocity (ignoring surface friction due to different landcover types - if this can be ignored) by choosing the steepest path.

The end application would be in determining the overland flow path from a spill site to a stream/river reach that provides the earliest possible time for entry of the spill into the stream/river.  There are different assumptions or starting datasets that could change the answer.  For instance, a 30m digital elevation model (DEM) cold prove useful in a forested or agricultural area, but a 1m or sub-1m DEM might be necessary for sorting out the true flow path (and therefore the true shortest path) in a highly urbanized area.  In either of these cases, however, once one has a DEM that he/she feels represents the terrain well enough for an analysis, which algorithm should they choose?

This query is borne out of trying to avoid calculating the distance from the spill to the stream/river reach using simple straight line geometry in a GIS tool, such as the proximity toolset in ArcMap.

Does anyone have any experience with this?

EDIT:  Based on the responses so far, it seems that using a routing algorithm that allows for dispersion, and then calculating the shortest path out of the possible paths predicted by the algorithm, might be the solution I want.  In the answer suggesting the use of TauDEM, the journal article that is linked to contains this image:

(image from Tesfa et al., 2011 - http://www.neng.usu.edu/cee/faculty/dtarb/Tesfa_EMS2011.pdf)

Examining the predicted flow path of the D-Inf algorithm, we see that several different paths are possible.  Using such an algorithm, I suppose on could calculate the shortest path, the longest path, and a range of other statistics on the predicted pathways.  If similar or near identical ground cover is found across the different pathways, then some combination of steepness and path length may dictate the quickest path?  Which might be different from the shortest path?  I'm starting to think I have to review my question and/or thought process - it may be necessary to incorporate equations/models governing the speed of overland flow along side the algorithms that examine possible flow pathways.  I need to keep reading and ingesting the article, but I think that the answer lies within the ideas presented in the article - so for now I will mark the suggesting answer as correct.


texts: Trigonometry and calculating the elevation angle from aerial photography

texts:
I was wondering whether anybody could advise me on a way to calculate the elevation angle from a camera to an object in an aerial photograph taken from a known elevation with no suitable reference objects.
Ultimately I would like to measure particular features within an image. I was thinking of a solution using trigonometry. Attached is a diagram illustrating the problem, as well as a sample of a picture to apply this to.
In the diagram, A represents the position of the camera. The Red triangle represents the focal view of the images. B is the centerpoint of the images. The Blue rectangles are objects I would like to measure.
To do it this way, I need to solve either angle E or F, however given that there is not tilt/angle/axis data stored within the photos EXIF metadata, I think the only way to do this would be to project the ground plane onto the image which I believe will allow me to calculate angle F. However I am unsure or unaware of whether there is software which may make this feasible?? As the goal would be for the solution to be applied to a series of photographs where the elevation angle could vary widely.
There aren't parallel lines, shadows, or objects of known size. However, topography shouldn't be an issue as the surface is flat (see image of ice).


texts: What's the difference between EPSG: 3043 and EPSG: 25831?

texts:
When I open a layer projected in ArcGIS with EPSG:25831, QGIS detect it as EPSG:3043.

Apparently everything is where it is supposed to be, and checking both EPSG:25831 and EPSG:3043 they appear to have the same bounds.

Is there any reason why QGIS deals with this info as 3043 or it is a bug? Is it possible that not changing it to 25831 can cause my data to be imprecise?


texts: What projection is being used in this dataset?

texts:
I'm not from US and so don't know much about the CRS/EPSG there. Anyway, when reviewing the FAA 5010 Airport dataset on https://www.faa.gov/airports/airport_safety/airportdata_5010/#5010 (direct download link = https://www.faa.gov/airports/airport_safety/airportdata_5010/menu/nfdcfacilitiesexport.cfm?Region=&District=&State=&County=&City=&Use=&Certification=) I found myself confused about what projection was being used. Can anyone help?

There is a "data dictionary" file that can be downloaded but CRS/EPSG is not included:

ARPLatitude = Airport reference point latitude (formatted).
ARPLatitudeS = Airport reference point latitude (seconds).
ARPLongitude = Airport reference point longitude (formatted).
ARPLongitudeS = Airport reference point longitude (seconds).

Here is a sample of the lat/long fields being used:

ARPLatitude |   ARPLatitudeS    | ARPLongitude |    ARPLongitudeS

51-53-00.8954N | 186780.8954N | 176-38-32.9277W | 635912.9277W

56-56-19.2870N | 204979.2870N | 154-10-57.2000W | 555057.2000W

60-54-49.7150N | 219289.7150N | 161-29-35.9850W | 581375.9850W

60-54-28.3130N | 219268.3130N | 161-26-06.2780W | 581166.2780W


texts:
The scale of my map is wrong, no matter what coordinate system I use in my project. The map and all objects show the correct coordinates but the scale of the map shows nonsense and the measurement tool gives wrong dimension.

The map should be like 1:10000 but qgis shows 1:1. Anyone an idea?


texts: How to prevent Mapbox GL JS FlyTo animation from being interrupted?

texts:
If a user clicks, swipe scrolls on mobile, or interacts with the map in any way, the flyTo animation is canceled. How can I stop that?

Or at the very least, jump to the destination?


texts:
I'm new to pyqgis and I'm not sure how can I find out if checkbox is checked and according to the result do different things in another method?

I know how to do it in theory, but if someone can show me your code I will be thankfull.


texts: Trying to loop through a directory and create a shapefile of a given threshold using arcpy.

texts:
Below is my code. I'm having issues with this will work sometimes then other times it will not work.

# Import dependencies:
import arcpy
import glob, os
arcpy.CheckOutExtension('Spatial')
arcpy.env.overwriteOutput = True

# Define variables:
filepath = r'\'
folder = []
threshold = -19

os.chdir(filepath)
for i in glob.glob('*.tif'):
    folder.append(i)

# Loop through the folder:
ras = arcpy.Raster(i)
body = arcpy.sa.Con(ras < threshold, 1)
outfile = r'' + i[17:35] + '.shp'
body_shp = arcpy.RasterToPolygon_conversion(body, outfile, "NO_SIMPLIFY")


texts: Is there a tool I could use to map Instagram photos that are using a specific tag? I don't necessarily want to show the image on a map, but rather a pin that someone could hover over to see the image. To add context, I'm developing a game and would like to include such a feature on the community website.


texts: Spatialite view appears 'filtered' in attribute table, won't load in canvas

texts:
I have a points table (all_trees) and polygons table (workareas) and a view (affected_trees_view) as follows:

SELECT * FROM all_trees, workareas
WHERE intersects(st_buffer(
                           all_trees.geometry,
                           CASE WHEN rpz_m IS NULL THEN 1 ELSE rpz_m END),
                 workareas.geometry)
GROUP BY all_trees.treeid

This returns 92 point features, which show up in DB manager/Qspatialite as expected.

However, when I try to load the view into QGIS I have two issues:

Points don't always visually appear on the canvas (if I load it as a spatial view using Qspatialite it works, sometimes... whether or not I manually insert the views_geometry_columns information).
Only see a single point in the attribute table even if the points visually appear on the canvas.

The feature count in Qspatialite, the layers panel, and attribute table status bar show the right number (92), but the attribute table shows "Filtered: 1" - see highlighted sections below

So I then created a table affected_trees from the view using

CREATE TABLE affected_trees AS SELECT * FROM affected_trees view;'

and that works perfectly fine - correct feature count, features appear, attribute table shows everything. (see second attribute table behind the first one below)

Why and how does the data get filtered like that?

It appears to me that the data itself is not exactly corrupted, it's just got something to do with how QGIS handles spatial views.

I've looked up whatever I could with my limited understanding of databases; the issue I have is somewhat different to this, for example. I've followed the latest answer from Can QGIS read Spatialite views? regarding updating views_geometry_columns - I've tried hand writing the view, I've tried using Qspatialite to generate the spatial view, I've used the following code:

INSERT INTO views_geometry_columns
(view_name, view_geometry, view_rowid, f_table_name, f_geometry_column,
read_only )

VALUES ('affected_trees_view', 'geometry', 'ROWID', 'all_trees', 'geometry', 1);

But I still encounter the same issues (particularly #2). I've posted a question about a similar issue with QGIS Virtual Layers here, though I have no idea if the two are connected.

I would like to be able to show the spatial view in QGIS so it shows points as soon as I change the dependent factors (workareas geometry, rpz_m value), rather than keep having to update the table with triggers.


texts: Large Legend doesn't appear on print task

texts:
I am having a problem displaying a large Legend (with more than 60 items) on the print task, it appears that the legend doesn't appear at all at the printed template.

It does print fine if there are few results (I tested a legend with less than 25 items and it works fine).

I am using the default Print widget of the webappbuilder tool with a little customization.

I am using ArcGIS JavaScript v3.18 and WebAppBuilder

Is there any solution to this ?

Here is a code sample of the point where the webappbuilder execute a print request.

 var template = new PrintTemplate();
 template.format = form.format;
 template.layout = form.layout;
 template.preserveScale = (form.preserveScale === 'true' || form.preserveScale === 'force');
 template.label = form.title;
 template.exportOptions = mapOnlyForm;
 template.layoutOptions = {
 authorText: hasAuthorText ? form.author : "",
 copyrightText: hasCopyrightText ? (form.copyright || this._getMapAttribution()) : "",
 legendLayers: legendLayers,//contains LegendLayer array that contains ids of the operational layer added to map.
 titleText: hasTitleText ? form.title : "",
                customTextElements: cteArray
            };
 this.printparams.template = template;
 this.printparams.extraParameters = { // come from source code of jsapi
                printFlag: true
            };
 var fileHandel = this.printTask.execute(this.printparams, function (res) {
          //res contains a link to the new printed template.
            });


texts: postgresql export database by command or pgAdmin III

texts:
I want to export my postgresql db with pgadmin.
I search all of its GUI, I also searched the web and I don't find some sueful information.
I know how to import the db into postgresql: create a db with `template0' and then:

 c:\program files\postgresql\9.3\bin> psql -h localhost -p 5432
                                -U postgres newDBname < D:\backup.sql

But I don't know how to export it.
please help me to export my db.


texts:
I'm trying to separate out records in a quasi-summarized table.

Currently, my features are geocoded address points. Each address point has one record associated with it with a field showing how many trees were planted at the property. Most of the time there will only be one tree per property but sometimes there are multiple trees planted at the property.

Ideally I'd like to separate/expand my records so that the number of attribute records for a particular address reflects the number of trees planted at that address.

For example, if property A has 3 tree plantings there will be 3 corresponding records in the attribute table (one-to-many).

If property B has 1 tree planting there will be 1 record in the attribute table.

Is there a tool or workflow that can automate this process??


texts: How to create a GeoJSON that works with D3?

texts:
I'm simply trying to convert a .shp file to a geoJSON format using:

ogr2ogr -f geoJSON output.json input.shp

After executing the command nothing seems to be wrong. Here's an excerpt from the output.json

    {
    "type": "FeatureCollection",

    "features": [
    { "type": "Feature",
    "properties": { "ID_0": 86, "ISO": "DEU", "NAME_0": "Germany", "ID_1": 1, "NAME_1": "Baden-Württemberg", "NL_NAME_1": null, "VARNAME_1": null, "TYPE_1": "Land", "ENGTYPE_1": "State" },
    "geometry": { "type": "MultiPolygon", "coordinates": [ [ [ [ 8.708400940398242, 47.715581894910606 ], [ 8.713716147005524, 47.701734382960055 ],
...

But when i try to use de JSON file in d3 (http://d3js.org/) to draw SVG polygons the output is just wrong. Since the shp files are displayed correctly in QGIS I think there has to be something wrong with the way i use ogr2ogr. The SVG I get is not completely wrong but there seems to be a detail I can't find. It seems like it has been turned upside down and somehow distorted into two separated parts.

Here's the javaScript i used to generate the svg:

//dimensions
var w = 2000;
var h = 2000;

var svg = d3.select("#chart").append("svg")
    .attr("width", w)
    .attr("height", h);

    d3.json(
    "http://localhost:8888/data/data.json",
    function (json) {

    var path = d3.geo.path();

    svg.append("g")
        .attr("class", "black")
        .selectAll("path")
        .data(json.features)
        .enter()
        .append("path")
        .attr("d", path);

Does anyone have an idea what went wrong here? I also tried to convert the shp-file using Qgis and myGeodata (http://converter.mygeodata.eu/vector). But neither of them work the way they should.

I'm very new to this whole cartography stuff. So I would be very happy to get some advice.

Thanks very much!


texts: How do I remove the terrain checkbox from Google Maps?

texts:
I don't mind having the map/satellite options, but I don't need the terrain option for my map. How do I disable just that? I don't see it documented well on

https://developers.google.com/maps/documentation/javascript/controls


texts: Debugging function takes exactly 3 arguments (1 given)?

texts:
Create a function called "numSchools" that counts the schools of a specific type.
The function should have three input parameters, (1) a string for the workspace, (2) a string for the shapefile name, and (3) a string for the facility type (e.g. "HIGH SCHOOL"), and one output parameter, (1) an integer for the number of schools of that facility type in the shapefile.

import arcpy

shapefile = "Schools.shp"
work = r"c:\Scripts\Lab 6 Data"
sTyp = "HIGH SCHOOL"

def numSchools(work, shapefile, sTyp):

    whereClause = "\"FACILITY\" = 'HIGH SCHOOL' " # where clause for high schools
    field = ['FACILITY']
    searchCurs = arcpy.SearchCursor(shapefile, field, whereClause)
    row = searchCurs.next()
    for row in searchCurs:
        # using getValue() to get the name of the high school
        value = row.getValue("NAME")

        high_schools = [row[0] for row in arcpy.SearchCursor(shapefile, field, whereClause)]
        count = arcpy.GetCount_management(high_schools)
        return count

numSchools("There are a total of: "),count

My error message states "  line 21, in
TypeError: numSchools() takes exactly 3 arguments (1 given)"
This error is for the final line of my code.

My code is only supposed to output the number of high schools.


texts: Joining a table to a polygon using Esri Insights for ArcGIS

texts:
I am trying to perform a many to one join using Esri Insights for ArcGIS. With in Insights I am trying to join a flow and pressure table to a polygon. Esri Insights does not allow for spatial joins or query tables so I am looking for a simple way to perform the join. The data format for the flow and pressure table is as follows Date, Time, Flow, Pressure and the field set up for the polygon is object id, shape, shape_Length, and shape_Area. The overall goal is to have the flow and pressure data appear in the attribute table of the polygon.

Is there a simple way I can format the data table or polygon to allow for the join ?


texts:
I have an extension which can be loaded in SQLite using the load_extension function. It would be useful for me if I could load it in QGIS to use in queries made via Database Manager.

The Python SQLite binding prohibits loading extensions by default. Could I do something from the Python console or project macro to load an extension?

I have QGIS 2.4.0 with spatialite 4.1.1, python-spatialite 3.0.1, python 2.7.8 and sqlite3 3.8.5. I have tried running the select load_extension('...') from DB Manager and it does not work (gives error "Authorization required"; because it is prohibited by default by python sqlite bindings). I don't know how to get the connection from the python console to do anything with it.


texts:
There was a very interesting feature in Map Book plugin for ArcGIS which was the possibility of not including in the resulting map book those pages where there's no feature from the layer or layers you want to map.

I've been reading QGIS documentation but I can't find an equivalent feature in Atlas tool. Is there any? Or maybe some trick to get this?


texts: Merging polygons to new shapefile with sum of total area using ArcMap

texts:
I am trying to group some polygons of a shapefile into a new one (with one record) but I want to maintain their total sum area. Some of them are overlapping. When I try Merge from Editor, the new total area is not the sum that I want.
How do I do this?


texts:
I have a rasters i want to publish it on Geoserver i do:

@echo off
set ext=.tif
set extout=.tif
FOR /f "tokens=*" %%G IN ('dir /b D:\OFP\*.tif') DO ("C:\Program Files\GDAL\gdalwarp" -r average -srcnodata 0  -s_srs "local srs proj" -t_srs EPSG:3857 D:\OFP\%%~nG%ext% D:\ofp_ave\%%~nG%extout%)

After publish as Image Mosaic, and got result:

There is white lines between rasters. I played with layer params:

But lines stand still.  When i open geotiffs in Qgis i do not see any lines.
How can i remove ut?


texts:
I want to multiply two raster layers. Values of the first layer range from 0 to 5 and values of the second layer range from 5 to 15. It is obvious that I should get some 0 values after multiplication, but MIN value of the output is 5. It means, that 0*5 is returning 5. I want it to return 0. I tried several things, for example:
 y3<-y2*y1
 y3<-overlay(y1,y2,fun=function(x,y){return(x*y)})

These return the same thing. I also tried to add arguments na.rm=TRUE to the function in second example, but nothing happened.


texts: gdal_rasterize using a "string" as src_datasource

texts:
I would like to use a string as src_datasource for gdal_rasterize program.

Is it possible to use something like a WKT definition to translate a single point to a raster?

Something like:

gdal_rasterize -b 1 -b 2 -b 3 -burn 255 -burn 0 -burn 0 -l mask "POINT(6 10)" work.tif


texts:
I'm trying to classify a set of buildings (grey) in a layer, according to construction year, in a number of different areas (orange).

In Orange – areas with a sum of buildings inside.

The construction year information is inside the area’s layer…
Ok, so I think this is a problem with data display, because all my data is organized like this:

So for Area D, there are 18 buildings built in 1940, and so on…

How can I classify this information to calculate maximum and show inside of each area what the most represented year is?


texts:
I am trying to use an update cursor on a shapefile to calculate a row from a dictionary and must be missing something very simple. below is my code

so my dictionary is called adict

adict = {'source': {0: u'SUN', 1: u'DAN'}, u'WORK_TYPE': {0: u'FRA', 1: u'UKB'}, u'API_NUM': {0: '422973', 1: '422873'}

And i am trying to calculate worktype in my shapefile.

with arcpy.da.UpdateCursor(myshapefile, ['API_NUM','WORK_TYPE']) as cursor:
    for row in cursor:
        if row[0] in adict['API_NUM'].values():
            row[1]=adict[row[0]]
            cursor.updateRow(row)
        else:
            print 'no'

i'm trying to set the row value on the shapefile to the corresponding dictionary value.

I am trying to calculate my work_type column in the shapefile with the work_type value from the dictionary based on the matching API_NUM, so for api = 422873, the worktype would be UKB.

If i were using geopandas, i would convert the dict to a dataframe and merge on api. Using arcpy right now i am iterating through both sets of data and calculating the work value from the corresponding dict data, but because there are so many records in both sets, it's taking a while using that method. so i thought a dictionary may be faster. I know i could convert to array, then to a table, but that just seems like a bunch of steps.


texts: Spelling of Kilometres on ArcGIS Pro scale bar?

texts:
Being in Australia I like to see the spelling of "Kilometres" used on my maps in preference to "Kilometers".

In ArcGIS Pro 1.0 I have not been able to find an option to change the text used for units and so my scale bar looks as below.

I have looked at the Element Properties for the Scale Bar and was hoping that setting Region and language to English (Australian) might do it but it did not.

Does anyone know how to set a spelling of "Kilometres" on an ArcGIS Pro scale bar?


texts: Generating tessellation in ArcGIS Desktop versions before 10.3

texts:
I'm using ArcGIS Desktop 10.2.2 and want to generate a tesselation grid but this tool is available in 10.3 or versions of ArcGIS.
Is there any way I can generate a hexagon grid for my study area using ArcGIS 10.2.2?


texts: Insert a position and check for geofence entry exit in postgres with postgis

texts:
Is it possible to insert a position, and any related geofence entry exits for that position in just one query in postgres (9.6) with postgis (2.3.2)?

For example given tables users, positions, fences and events:

user_id | name
1       | bob
2       | jane

position_id | user_id | lat        | lon
1           | 2       | 33.1261983 | -117.3138087

fence_id | geometry
1        | 0101000021E610000096F92563F71F63402664979471963EC1
2        | 0101000022E610000096F92563F71F63402664979471963EC2

event_id | type  | position_id | fence_id
1        | entry | 1           | 1
2        | entry | 1           | 2
3        | exit  | 1           | 3

I'm currently doing it in multiple steps:

Insert the position:

insert into positions (user_id, lat, lon) values ($1, $2, $3)
returning id as new_position_id

Get the fences for the new position and the previous position, with $2 param being the value returned from the insert above

with current_and_previous_position as (
    select id from positions where user_id = $1 and id < $2
    union all (select $2)
)
select * from
positions
join fences on ST_Contains(geometry, ST_SetSRID(ST_Point(lon, lat), 4326))
where positions.id in (select id from current_and_previous_position)

Then the next step is some Node.js code that takes the fences list and calculates if there are any entries or exits:

var entries = getEntriesFor(user, position, positionFences);
var exits = getExitsFor(user, position, positionFences);

And finally, I save the entries and exits (where $1 is the entry + exit lists concatenated together):

insert into events (type, position_id, fence_id)
select * from json_populate_recordset(null::events, $1)

This works, but takes too many individual steps, and 3 database operations for each and every position received.

Is it possible to do this in one single DB operation that inserts both the position, and the geofence entries and exits?


texts:
I have two rasters and am using QGIS 3.6

I want to know which areas in raster 1 that are "no data" are not "no data" in raster 2.

I also want to know what areas in raster 2 that are "no data" are not "no data" in raster 1.

I assume that the QGIS raster calculator is the tool to use. I dont think there is a way to use conditional statements and I am unsure how to enter No data into it.


texts: Why are Galileo satellites inaccessible in USA?

texts:
I see that some (or all?) Android devices that support Galileo global positioning system are deliberately programmed to self-disable their Galileo support when used in USA or US territories.

This applies to Sony Xperia XZ Premium, for one example, whose hardware supports GPS, GLONASS, BDS and Galileo. According to the official response, it is not a matter of region-specific firmware, but rather a purposely programmed real-time behavior: any XZ Premium device, regardless of its region of origin, will see Galileo satellites outside of USA, but inside USA it will deliberately pretend that they don't exist.

What is the reason for this behavior?

Is this some sort of company-specific patent or legal issue?

Or is it some kind of US or European regulation?


texts: Converting Points shapefile to 3D polylines in QGIS (without Grass)

texts:
I'm using QGIS 3.10 and am trying to come up with a way to create a shapefile of polylines with Z attribute attached, starting from a points shapefiles.

Both Saga/Convert Points to Lines and QGIS/Points to Path don't keep the Z values on the resulting polylines' vertices.

I've look up for hours online and I've only seen GRASS v.drape as a solution but GRASS hasn't been set up on my computer because of IT internal policies.

Has anyone any ideas about this?


texts: Finding centroid of features within polygon using ArcGIS Pro

texts:
I have some polygon features and within those polygons are some coordinates for some other features.
How do I find the centroid of those coordinate points within the polygon (NOT the geometric centroid of the TAZ)?


texts: What are properties of tile geometry of Sentinel data?

texts:
I have downloaded sample tile from

http://sentinel-s2-l1c.s3-website.eu-central-1.amazonaws.com/#tiles/10/S/DG/2015/12/7/0/

and found tileInfo.json file with the following information

"tileGeometry" : {
    "type" : "Polygon",
    "crs" : {
      "type" : "name",
      "properties" : {
        "name" : "urn:ogc:def:crs:EPSG:8.8.1:32610"
      }
    },
    "coordinates" : [ [ [ 399960.0, 4200000.0 ], [ 509760.0, 4200000.0 ], [ 509760.0, 4090200.0 ], [ 399960.0, 4090200.0 ], [ 399960.0, 4200000.0 ] ] ]
  },

It appears, that tile is strictly rectangular in some coordinates, because corners contain duplicate values. Is this property constraint, i.e. true for all tiles?

What units are used for these coordinates? I tried to interpret them as lat-lon multiplied by some factor, but probably failed (was unable to identify any ground details with compre to google maps).

Also I can interpret it as meters, because they have 6 digits and meters should have 5.

UPDATE

I have plotted coordinates

Both height and width are 109800 of units, so it can't be meters, which would sum up to 100000. But it coincides with the number of pixels in tile.


texts: OpenLayers: GeoJSON Vector-Layer refresh / redraw every second

texts:
I do not get to move the Vector Layer in my OpenLayers script. The layer should be renewed every second. I have already tried a lot, like refresh() and redraw() or OpenLayers.Strategy.Refresh, but nothing has worked for me.

The GeoJSON-"polygon" should be re-layered every second.

Here is a JSFiddle:
https://jsfiddle.net/ejts8foc/4/