texts: What logical function should i use to do a union of two raster files in GRASS?

texts:
The files are :
DEM and roads. For the following i did euclid distance and then reclass. Now i need to know with what logical function to combine them

                   elevation   rank
               DEM < 300        10
               DEM > 300        0

                     distance    rank
proximity to roads < 2000        10
proximity to roads > 2000        5


texts: How to change page size in QGIS Composer

texts: I've made a few maps in QGIS Composer and set them up initially as 8.5x11 page size.  I want to begin with that layout and apply it to 11x17 page size.  However when I change this in the Page Layout, only the page dimensions change.  The symbol sizes should have changed as well but they look the same as if they were printed on an 8.5x11.  The fonts also remain the same size in relation to the composer page.  When the page is saved to PDF, Acrobat does not recognize that the page is 11x17 and still shows it as 8.5x11.  Am I missing something?


texts: How do I specify the QGIS specific Python version to add Python packages to that version?

texts: On my MacBook (running Mojave), I have preexisting Python versions (Anaconda and the Python 2.7 that comes with the Mac). I installed QGIS version 3.10.5-A Coruña. To use a QGIS plug-in (Enmap-Box), I need to add some Python packages. My Mac terminal window doesn't recognize the QGIS Python. What do I need to do to switch to the QGIS Python so that I can add packages there?


texts: Listing/Referencing Specific Bookmarks in Python script tool?

texts:
I'm trying to create three map layouts and export them into PDFs using ModelBuilder in ArcGIS Pro. However, when I try to reference the three bookmarks (AAA, Amityville, and Astoria) to be used for the layouts, I get an error.

Here is a section of my code from a custom script tool that I am having trouble with:

# setting layout to bookmark
aprx = arcpy.mp.ArcGISProject("Current")

# add name of layout
lyt = aprx.listLayouts("Layout_King")[0]
mf = lyt.listElements("MAPFRAME_ELEMENT")[0]

# add name of bookmark
bkmks = mf.map.listBookmarks('AAA','Amityville','Astoria')
bkmks.sort(key=lambda x: x.name, reverse=True)
for bkmk in bkmks:
  mf.zoomToBookmark(bkmk)
  legend.autoAdd = True
  time.sleep(5)
  lyt.exportToPDF(r"C:\arcGIS_Shared\Exports" + "\\" + bkmk.name + ".pdf")

The error I'm getting is:

Traceback (most recent call last):
  File "C:\arcGIS_Shared\PythonTest_2\Layout_Bookmark_Export_Script.py", line 42, in <module>
    bkmks = mf.map.listBookmarks('AAA','Amityville','Astoria')
TypeError: listBookmarks() takes from 1 to 2 positional arguments but 4 were given
Failed to execute (Layout Bookmark Export)


texts: Long loading time while opening layer preferences

texts:
Recently I have changed my HDD to a SSD type, then reinstalled QGis and loaded the project I was working on. Everything works as good as it should, with only one little problem:
The loading time of layer preferences (double clicking one specific layer) takes 15 to 20 seconds, while the whole program freezes. Once the preferences window opens everything works fine again.
Does anyone had a similiar issue or has an idea why this loading time might occur?
Thanks alot in advance!


texts: Clear Environment Settings in scratchWorkspace

texts:
I have several mxd files with a bad default geo-processing environment. The current and scratch workspaces are set to a location that no longer exists. I'm trying to use python to clear those settings and reset to ArcGIS default, which is blank, so that the user's environment or default geodatabase will be used. My python script successfully clears the info, but does not save it. How can I fix my script to save the default environment for each mxd? I'm missing something.

import arcpy, os
Workspace = r"..."
arcpy.env.workspace = Workspace
mxdList = arcpy.ListFiles("*.mxd")
print(mxdList)
for mxd in mxdList:
    filePath = os.path.join(Workspace, mxd)
    mapedit = arcpy.mapping.MapDocument(filePath)
    arcpy.ResetEnvironments()
    arcpy.ClearEnvironment("workspace")
    arcpy.ClearEnvironment("scratchWorkspace")
    print(arcpy.env.workspace)
    print(arcpy.env.scratchWorkspace)
    mapedit.save()
    del mxd


texts: Adding additional points to GCP table after an image has been georeferenced

texts:
I have finished georeferencing an image with no attached coordinate system. I was unsatisfied as the image was not accurately overlayed on the map. The dx/dy & Residual error values are not significantly high indicating that there was no human error in the georeferencing process.

I am trying to re-open the GCP table to add more points, however, I have been unable to find a way to do this. The only solution I can think of is to re-do  the  georeferencing. Is there a way to re-open the GCP table to add more points?

Version: QGIS 3.6
Operating System: Windows 10


texts:
I have a two shapefiles, one is a DMA level and one of US states level.

How do I create a shapefile containing DMA of a certain state, for instance, Texas from a larger shapefile.

I tried running this query, but with no success.

SELECT geom FROM usdma
WHERE ST_Overlaps(usdma.geom, usstates.geom)
WHERE usstates.name = 'TX';


texts: Clipping MrSID raster image, without decompressing clipped raster, in ArcGIS Desktop?

texts:
I am using ArcView ArcMap10 and when I clip a raster using tools/data management tools/raster/clip. The Clipped file ends up much larger than the original file.  Specifically I am trying to clip a MrSID file that is 2.6gb, the clipped file ends up being 54.gb.  I have tried playing with the environment setting and compression rates trying to convert to a JP2 with no luck.

I am fairly new to GIS.


texts:
I use postgresql 8.4 to route a river network, and I want to use psycopg2 to loop through all data points in my river network.

#set up python and postgresql connection
import psycopg2

query = """
    select *
    from driving_distance ($$
        select
            gid as id,
            start_id::int4 as source,
            end_id::int4 as target,
            shape_leng::double precision as cost
        from network
        $$, %s, %s, %s, %s
    )
;"""

conn = psycopg2.connect("dbname = 'routing_template' user = 'postgres' host = 'localhost' password = '****'")
cur = conn.cursor()
while True:
    i = 1
    if i <= 2:
        cur.execute(query, (i, 1000000, False, False))
        i = i + 1
    else:
        break
rs = cur.fetchall()
conn.close()
print rs

The code above costs a lot of time to run even though I have set the maximum iterator i equals to 2, and the output is an error message contains garbage,

I am thinking that if postgresql can accept only one result at one time, so I tried to put this line in my loop,

rs(i) = cur.fetchall()

and the error message said that this line has bugs,

I know that I can't write code like rs(i), but I don't know the replacement to validate my assumption.

So should I save one result to a file first then use the next iterator to run the loop, and again and again?

I am working with postgresql 8.4, python 2.7.6 under Windows 8.1 x64.

Update#1

I can do loop using the code below, and the result follows,

import sys
import psycopg2

query = """
    select *
    from driving_distance ($$
        select
            gid as id,
            start_id::int4 as source,
            end_id::int4 as target,
            shape_leng::double precision as cost
        from network
        $$, %s, %s, %s, %s
    )
;"""

conn = psycopg2.connect("dbname = 'routing_template' user = 'postgres' host = 'localhost' password = '****'")
cur = conn.cursor()

rs = []
i = 1
while i <= k:
    cur.execute(query, (i, 1000000, False, False))
    rs.append(cur.fetchall())
    i = i + 1

for record in rs:
    print record

conn.close()

[(1, 2, 0.0), (2, 2, 4729.33082850235), (3, 19, 4874.27571718902), (4, 3, 7397.215962901), (5, 4,
6640.31749097187), (6, 7, 10285.3869655786), (7, 7, 14376.1087618696), (8, 5, 15053.164236979), (9, 10, 16243.5973710466), (10, 8, 19307.3024368889), (11, 9, 21654.8669532788), (12, 11, 23522.6224229233), (13, 18, 29706.6964721152), (14, 21, 24034.6792693279), (15, 18, 25408.306370489), (16, 20, 34204.1769580924), (17, 11, 26465.8348728118), (18, 20, 38596.7313209197), (19, 13, 35184.9925532175), (20, 16, 36530.059646027), (21, 15, 35789.4069722436), (22, 15, 38168.1750567026)]
[(1, 2, 4729.33082850235), (2, 2, 0.0), (3, 19, 144.944888686669), (4, 3, 2667.88513439865), (5, 4, 1910.98666246952), (6, 7, 5556.05613707624), (7, 7, 9646.77793336723), (8, 5, 10323.8334084767), (9, 10, 11514.2665425442), (10, 8, 14577.9716083866), (11, 9, 16925.5361247765), (12, 11, 18793.2915944209), (13, 18, 24977.3656436129), (14, 21, 19305.3484408255), (15, 18, 20678.9755419867), (16, 20, 29474.8461295901), (17, 11, 21736.5040443094), (18, 20, 33867.4004924174), (19, 13, 30455.6617247151), (20, 16, 31800.7288175247), (21, 15, 31060.0761437413), (22, 15, 33438.8442282003)]

but if I want to get this look of output,

(1, 2, 7397.215962901)
(2, 2, 2667.88513439865)
(3, 19, 2522.94024571198)
(4, 3, 0.0)
(5, 4, 4288.98201949483)
(6, 7, 7934.05149410155)
(7, 7, 12024.7732903925)
(8, 5, 12701.828765502)
(9, 10, 13892.2618995696)
(10, 8, 16955.9669654119)
(11, 9, 19303.5314818018)
(12, 11, 21171.2869514462)
(13, 18, 27355.3610006382)
(14, 21, 21683.3437978508)
(15, 18, 23056.970899012)
(16, 20, 31852.8414866154)
(17, 11, 24114.4994013347)
(18, 20, 36245.3958494427)
(19, 13, 32833.6570817404)
(20, 16, 34178.72417455)
(21, 15, 33438.0715007666)
(22, 15, 35816.8395852256)

What should I make a little change in the code？


texts: How to add features where there are no other features in Openlayers 4.x

texts:
I would like to add features to a map view where there are no other features within a certain radius.

I am guessing I have to use some sort of R-tree or Quad-tree search of the view extent to find points that would fit that requirement. Or perhaps some packing algorithm? Is there some simpler way using the API?


texts: Watching for scale/zoom changes in ArcGIS API for JavaScript 4?

texts:
I am trying to add an event listener for when the scale of the map changes, that is scale or zoom level changes in order to enable or disable feature layer label text. When the user zooms in to a larger scale I want the labels to be visible and when the scale is lowered the labels should go away. After defining the MapView I tried the following

 view.on("zoom", function(){alert();});


texts: Changing file names inside the shapefile.zip file obtained from a WFS

texts:
I am using the following format to download a shapefile from a GeoServer WFS.

http://some.domain.com/geoserver/WFC/ows?service=WFS&version=1.0.0&request=GetFeature&typeName=WFC:someLayerName&maxfeatures=5000&outputformat=SHAPE-ZIP

I receive a shapefile named WFC:someLayerName.zip and the different files inside that zip folder are named as someLayerNamePolygon.shp or someLayerNamePolygon.cst etc.
I was wondering if there is a way to change the name of files INSIDE the zip folder? I am aware that the name of the actual shapefile.zip can be changed but that does not affect the name of the files inside the zip folder.


texts:
I'm still new in the world of GIS.

I'm currently doing some zonation with "test points."  Each has its own coordinates and "P" value. I want to make contours of these data, to divide the whole area into three zones (Yellow zone if P <= 5, Orange zone if 5 < P <= 15, and Red zone if P > 15).

Please guide me step to step on how to make this with ArcGIS 10 (I don't have a license for 3D Analyst).


texts: Problem executing Processing Vector Grid algorithm in a standalone script

texts:
I have the following piece of code which I execute within QGIS 2.8.6

import processing
import ogr
restrLayer = QgsVectorLayer("/home/usr/Desktop/trial/EligibleAreas_polygons_filtered.shp", "EligibleAreas_polygons_filtered", "ogr")
if not restrLayer.isValid():
    print "Restrictions Layer Failed To Load"
QgsMapLayerRegistry.instance().addMapLayer(restrLayer)
cellsize = 1000
xmin = (restrLayer.extent().xMinimum())
xmax = (restrLayer.extent().xMaximum())
ymin = (restrLayer.extent().yMinimum())
ymax = (restrLayer.extent().yMaximum())

extent = str(xmin)+ ',' + str(xmax)+ ',' +str(ymin)+ ',' +str(ymax)
grid1 = "/home/usr/Desktop/trial/grid1.shp"
processing.runalg('qgis:vectorgrid', extent, cellsize, cellsize, 0, grid1)

and it works fine. When I create a standalone script and try to do the exact same thing outside of QGIS I have the following code:

import qgis
from osgeo import gdal
from osgeo import ogr

from gdalconst import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *
from qgis.core import *
from qgis.analysis import *
from qgis.gui import *

import sys
import glob
import re
import os, fnmatch
import numpy as np
import numpy.ma as ma
import math

qgisprefix = '/usr'

# configure paths for QGIS
os.environ['PATH'] = qgisprefix+'/bin'
os.environ['LD_LIBRARY_PATH'] = qgisprefix+'/lib'
sys.path.insert(0, qgisprefix+'/share/qgis/python')
sys.path.insert(1, qgisprefix+'/share/qgis/python/plugins')

# disable QGIS debug messages
os.environ['QGIS_DEBUG'] = '-1'

# configure QGIS paths
QgsApplication.setPrefixPath(qgisprefix, True)

# initalise QGIS (load providers)
QgsApplication.initQgis()

# create a new application object
# without this importing processing causes the following error:
# QPixmap: Must construct a QApplication before a QPaintDevice
app = QgsApplication([], True)

from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *
print Processing.getAlgorithm("qgis:vectorgrid")

#add the restrictions layer
restrLayer = QgsVectorLayer("/home/usr/Desktop/trial/EligibleAreas_polygons_filtered.shp", "EligibleAreas_polygons_filtered", "ogr")
if not restrLayer.isValid():
   print "Restrictions Layer Failed To Load"
QgsMapLayerRegistry.instance().addMapLayer(restrLayer)

cellsize = 1000

xmin = (restrLayer.extent().xMinimum())
xmax = (restrLayer.extent().xMaximum())
ymin = (restrLayer.extent().yMinimum())
ymax = (restrLayer.extent().yMaximum())

extent = str(xmin)+ ',' + str(xmax)+ ',' +str(ymin)+ ',' +str(ymax)

grid1 = "/home/usr/Desktop/trial/grid1.shp"
processing.runalg('qgis:vectorgrid', extent, cellsize, cellsize, 0, grid1)

I get no output, so no grid1.shp. From print Processing.getAlgorithm("qgis:vectorgrid")  I get:
ALGORITHM: Vector grid
    EXTENT <ParameterExtent>
    STEP_X <ParameterNumber>
    STEP_Y <ParameterNumber>
    TYPE <ParameterSelection>
    OUTPUT <OutputVector>

so the needed algorithm from processing is found but for some reason the output layer is not created. I also checked the extent I get from the standalone script and it is the same as what I have in the python console (so the layer is loaded correctly). What could be the reason why I am not getting an output?


texts:
In QGIS I wanted to move points to the nearest lines. I checked the solutions here, but there are using tools which I am not familiar to yet (PostGIS or GRASS GIS). I figured out my own way with using QGIS but it is giving me weird results.
My problem is that moved points are not exactly on the lines, there is 0.0000 length line beetwen them.
Working process: I connected the lines with points using v.net.connect, then I extracted only the lines beetween points and orignal lines layer (let's call them join). I created start and ending points on the join layer using small python code:
from qgis.core import *
import qgis.utils

layers = iface.legendInterface().layers()

points = []

for i in layers[2].getFeatures():
    geom = i.geometry()
    a = geom.asPolyline()
    points.append(a)

vl = QgsVectorLayer("Point", "beg_end", "memory")
pr = vl.dataProvider()

vl.startEditing()

for i in range(len(points)-1):
    fet = QgsFeature()
    fet2 = QgsFeature()
    fet.setGeometry(QgsGeometry.fromPoint(QgsPoint(points[i][0])))
    fet2.setGeometry(QgsGeometry.fromPoint(QgsPoint(points[i][len(points[i])-1])))
    pr.addFeatures([fet])
    pr.addFeatures([fet2])

vl.commitChanges()
vl.updateExtents()

QgsMapLayerRegistry.instance().addMapLayer(vl)

In the end I saved the end points of the join which are supposed to lie on my line layer (and they look like they are there). But if I use connect tool from networks plugin it turns out that it is creating 0.0000 length line between moved points and original lines to join them.
Why the points are not on the lines and what is the reason of this small differences? Is it possible to snap them to the lines and how?


texts: What is the method for finding the percentage of common cells among 2 or more DTMs for the same area?

texts: I have a study area and 6 different DTMs for that area from different sources and I would like to know the percentage of pixels with the same value. I suppose I will have to extract some statistic and make some calculations in Excel, but I don't know what to extract and how to do it. If somebody knows some methods or something similar please give some indications.


texts:
How do Pix4D and Agisoft Photoscan mechanically differ in producing an orthomosaic, given the same set of geotagged images? Are they any resources that discussed this? I'm also assuming full capabilities for each software.

As I understand it, in a general sense they both consist of calculating camera positions, creating a point cloud, a digital elevation model (DEM), and then an orthomosaic. However, I didn't know specifically if they use different steps or mathematics to get this final product.


texts:
Recently I am working on a project which is based on making map for state-wise health care indicators. And I am getting the data (shapefile) from here: link
I just downloaded the shapefile which are present in that site and then I tried to add that shapefile in QGIS software but I can't and it shows :
Invalid Data Source: C:\Users\sasua\Desktop\India_District_Boundary.shp
is not a valid or recognized data source

I am pretty sure that file is shapefile but it doesn't work. What do I have to do to solve this problem?


texts: Reprojecting a population raster from LAEA to lonlat

texts:
I have a population density raster from the EEA which is based on the CORINE land cover 2000 data.

http://www.eea.europa.eu/data-and-maps/data/population-density-disaggregated-with-corine-land-cover-2000-2

The population density is given in inhabitants/km^2, the resolution is 100mx100m.

My goal is:

extract a country from the data by means of a shapefile from gadm.org
convert it to a population COUNT raster
reproject it from laea/GRS80 to longlat/wgs84
aggregate/resample it to match climate data resolution (0.11°)

I do the first step with QGIS, I load the whole .tif, then select the shapefile as mask.
Then I use the rastercalculator to divide by 100 to get the population count per cell.

First weird thing: R (using raster/rgdal packages) tells me that the minimum value is -1 (or -0.01) for the clipped (or clipped count) file, which is strange for a population raster.

The reprojection part doesn't work at all, as I understood, I'd use ngb for categorical data which I don't have. Bilinear resampling gives me a lot of negative values though.

The same problem occurs once the data is in the correct coordinate system, but not in the right resolution.

I first used the aggregate function in R to get the data close to the desired resolution, then used the resample function

Which is the correct reprojection method? How can I preserve the integrity of the population count?

Original data is in
+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

and should be transformed to
+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0


texts: Can you save a project using an ArcGIS Online Web App as a user without an account?

texts:
Once I make my web app publicly accessible, is there a way for users to save a project they're working on?

For example, an agricultural planning adviser is making a map of a farm they're visiting. They've drawn some polygons of the farm, made buffers around streams etc. - Can they save their work at this point and return to it at a later date?


texts: How can I constrain my plot to the extent of the raster image?

texts:
I am currently trying to create a multi-panel plot to display six raster images (3,2) using par() and plot(raster()). When I plot the rasters though, R is displaying the Y bounds far beyond what is necessary, resulting in lots of wasted space. How can I adjust this? When, I run extent() on the raster images themselves, the extent is correct.
a<-raster("file.tif")
...and so on to f.

par(mfrow=c(2,3))
plot(a)
plot(b)
plot(c)
plot(d)
plot(e)
plot(f)


texts: Points disappear when zooming past certain scale in ArcMap?

texts:
My points are disappearing when I zoom in past 300,000 in ArcMap.  I tried to repair geometry, but it doesn't fix anything?

I also checked the scale extents and those are turned off.  I can export the points to a new shapefile/feature class, but I'd like to not have to do that.

I'm using 10.1 on Windows 7 x64.


texts: Using Calculate Field for If statement over many fields?

texts:
I'm struggling with the python syntax in the Calculate Field so that, in the table above:

if b, c, d, or e are NOT NULL populate the Summary field with "Verified"
else populate Summary field with "In Progress"

So if any of the B, C, D or E columns contain NULL the Summary field will read "In Progress".

I believe this should be something along the lines of:

def Reclass !SUMMARY!:
    if [ !B! , !C! , !D! , !E!]  <> NULL
        !SUMMARY! = "Verified: "
    else
        !SUMMARY! = "In Progress: "


texts: Extruding mesh along parametric curve and applying earth curvature?

texts:
I'm working on a railway generation project, I must generate 3D meshes for the track (rail, ballast, ...). Track origin is defined geographically (latitude/longitude/azimuth/elevation), then the trajectory is defined using curvature (lateral and vertical) values. These curvature values permits to draw the track with three sort of parametric curve : straight lines, circle arc and clothoids (https://en.wikipedia.org/wiki/Euler_spiral).

For now I have a solution working correctly for small tracks, which consists of taking an origin point on the earth (near the track origin) and taking a tangent plane at this location (https://en.wikipedia.org/wiki/Local_tangent_plane_coordinates). Then all generation is done on this plane.

But, as one can imagine, on long tracks (several hundreds of kilometers), the meshes don't follow earth curvature. It is thus impossible to match real life datas.

What could be the best solution to address this earth curvature constraint?

Should we always work in spherical coordinate (and thus convert euclidean line/circle/clothoid parametric equation into spherical ones)?

Or could we apply a post transformation to remap meshes on earth?


texts: I am trying to use QGIS to assess what areas will flood if the river floods a specific height.  I have a DEM of the area as well as a vector file of the lines of the river of interest.  I cannot get a raster file of the river height.  Is there any way to do this?


texts: Creating contour lines ERROR 1 CPLRealloc() out of memory allocationg 416 bytes

texts:
When i try to create contour lines from raster (tiff file) in QGIS 2.6.0 Brighton i get the following message:

I thought it could be irregular size of raster that causes error so i clipped it in regular size but it did not help. I also tried to convert raster but it did not help either. When i take some other raster i create contour lines with no problems.


texts: Attribute join on the fly for visualization with ArcGIS JS API

texts:
Lets say I have 2 layers - one that contains shape information about a given group of entities, say countries, and one that contains some data attributes about the same group of entities.  Both layers have an identifier attribute that can be used to join the data, say, the country name, or ISO code.  To make this more concrete, let's take some examples.  On ArcGIS Online, we have the Countries reference layer, which has all the shape/polygon information, along with country name and ISO codes.  Renewable Energy as % of Consumption is a good example of the second type of layer - it has identifier attributes for county name and ISO code, as well as data points for energy consumption, but its only has spatial attributes about the country's latlng centroid, not any polygon data.
I want to display the countries as polygons, but colored according to the data in the second dataset.  I want to do this in an arcgis-js-api app. If all the data was part of one dataset, this would be easy. I could just do something like
const renderer = {
  type: 'simple',
  symbol: defaultSymbol,
  visualVariables: [
    {
      type: 'color',
      field: 'value_2017',
      stops: [
        {
          value: 0,
          color: '#ebe6df',
          label: '0',
        },
        {
          value: 500,
          color: '#7bccc4',
          label: '500',
        },
        {
          value: 11000,
          color: '#436480',
          label: '11,000',
        },
      ],
    },
  ],
};

export const energyByCountry = new FeatureLayer({
    url: countriesReferenceUrl,
    outFields: ['*'],
    renderer
});

Of course the issue with this is that the value_2017 attribute is not on the countriesReferenceUrl layer.  How can I tell the renderer to use the attribute value of the other layer for the feature that it matches to according to country name / ISO code / some ID value.
I had a similar issue last month, where I asked the question Forking and adding to FeatureLayer using ArcGIS Online.  This scenario was slightly different, as the second data source was not a featurelayer, but rather just an excel table with 0 geographic information, but it did have an ID column (for state names and abbreviations) that I could use to match the data with a polygon layer.  The solution here was to create a new layer - starting from the layer which contains polygons, I could do a data add-join in ArcGIS Pro, tacking on additional attributes to the existing shapefile layer.  I feel like in the case where both layers already exist as featurelayers, and there are matching ID fields, this seems like it would be overkill.
Is there a way in the ArcGIS API for JavaScript to do this on the fly or do I have to create a new layer which joins the data between the two existing layers?


texts: Does the MODIStsp package computes the QA's and eventually corrects the bad quality pixels?

texts:
I am using R's MODIStsp package for downloading MOD11A2 time series products.

While executing the command it was written that QA's were computed. Several files were created, LST_Day_1km, QC_Day, QAday_LST_err, QAday_mand, QAday_qual.

My question is, in which folder the 'corrected' images are? Or the package simply computes the QA's and I have to use the other package in order to produce images which contain pixels with the best quality?


texts: Transitioning from programmer to mapper?

texts:
I have now been rejected on 5 jobs asking for mappers on the pure basis that I am a programmer.  I don't get interviews and when I ask I get response of "You are a programmer and we wanted a mapper with programming experience".

I have my basics, I can use ArcMap, configure ArcSDE, understand coordinate systems.  I have experience writing ESRI add-in code in ArcMap and direct code accessing the server without ArcMap.  I understand how the spatial type is implemented and know how to create a matching UDT in accordance to OGC.  I am familiar with the basic spatial operations.

I know that I am missing some things (possibly a lot of things), but nobody wants to tell me what they are.  So please tell me.  Please assume some generality, I want to know what most programmers lack when compared to mappers, just listed my case to prove motive.

Note:  Open source GIS is coming next.


texts: Outputing shape_Area value from polygon feature to output table in ModelBuilder?

texts:
In ModelBuilder, I'm using the "Iterate Datasets" tool to iterate though a dataset of polygon features. Each are simple, single, polygons so just the one row in the attribute table. The polygons were generated with a shape_Length and shape_Area column in their attribute table.

I would like to output the area of each polygon feature to a summary table (cvs or similar) containing values for all iterations.

My current model is stuck right at the start with the Iterate Datasets tool which uses a Feature Dataset of polygon features as input. Previously I've just manually copied the shape_Area value from the features into an excel sheet which works but is slow.


texts:
I have a layer with several polygons and a layer with line above these polygons. Its like a map and a road network on it.
I have to find out the length of line in each polygon.
I used the intersect tool to spit the lines. However, I have three roads(from same layer) and when i intersect it. Instead of getting three different attribute, I get three lines but one single line attribute.
Could someone tell me how to measure the length of each line.


texts: QGIS 2.4 crashes when symbolizing QgsComposerShape in Python with ScriptRunner

texts:
I'm attempting to use ScriptRunner to output a jpeg of a red rectangle using QGIS Composer to understand the QgsComposerShape.  Everything works except QGIS crashes when the symbol is applied.  The output image is correctly produced but QGIS crashes.  If I remove the symbol portion of the code it draws the rectangle outline and doesn't crash.  If I run the same code line by line in the QGIS Python console everything works without crashing.

I just don't understand why using the ScriptRunner plug-in results in a crash.  Looking around GIS-SE I can't find this exact problem but there seems to be other issues with plug-ins and symbols.  Is there some workaround for this issue?

Here is the ScriptRunner script.  The layer isn't rendered but you need a map renderer in order to use the shape object:

from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *

def run_script(iface):
    """Print a simple map with a rectangle"""
    lyr = QgsVectorLayer("/Users/joellawhead/qgis_data/ms/Mississippi.shp", "Mississippi", "ogr")
    reg = QgsMapLayerRegistry.instance()
    reg.addMapLayer(lyr)
    c = QgsComposition(iface.mapCanvas().mapRenderer())
    c.setPlotStyle(QgsComposition.Print)
    c.setPaperSize(215.9, 279.4)
    w, h = c.paperWidth(), c.paperHeight()
    x = 0
    y = ((h * .2) * -1) / 2
    composerMap = QgsComposerMap(c,x,y,w,h)
    c.addItem(composerMap)
    dpi = c.printResolution()
    c.setPrintResolution(dpi)
    dpmm = dpi / 25.4
    width = int(dpmm * c.paperWidth())
    height = int(dpmm * c.paperHeight())

    red = {'color':'255,0,0,255','color_border':'0,0,255,255'}
    redsym = QgsFillSymbolV2.createSimple(red)
    shape1 = QgsComposerShape(10,50,10,25,c)
    shape1.setShapeType(1)
    shape1.setUseSymbolV2(True)
    shape1.setShapeStyleSymbol(redsym)
    c.addItem(shape1)

    # create output image and initialize it
    image = QImage(QSize(width, height), QImage.Format_ARGB32)
    image.setDotsPerMeterX(dpmm * 1000)
    image.setDotsPerMeterY(dpmm * 1000)
    image.fill(0)

    # render the composition
    imagePainter = QPainter(image)
    sourceArea = QRectF(0, 0, c.paperWidth(), c.paperHeight())
    targetArea = QRectF(0, 0, width, height)
    c.render(imagePainter, targetArea, sourceArea)
    imagePainter.end()

    image.save("/Users/joellawhead/qgis_data/map.jpg", "jpg")


texts: Set GPLinearUnit type parameter default based on input spatial reference (Buffer Example)

texts:
Within a .pyt python toolbox I'd like to have a parameter of data type GPLinearUnit get a default value based on an input Feature Class (FC). The standard Buffer (Analysis) tool works as a great example of what I'm trying to imitate:

The closest I've been able to come is to set units in updateParameters based on the FC spatialReference.linearUnitName, but linearUnitName and GPLinearUnit units aren't the same and there are 50+ linear units. Anyone know a better way to do this or a reference/code somewhere I could at least pull the linearUnitName to GPLinearUnit from?

def getParameterInfo(self):
    inFC = arcpy.Parameter(
        displayname = "Input FC",
        name = "inputFC",
        datatype = "DEFeatureClass",
        parameterType = "Required",
        direction = "input")
    units = arcpy.Parameter(
        displayName = "Distance & Units",
        name="Units",
        datatype="GPLinearUnit",
        parameterType="Required",
        direction = "Input")
    params = [inFC, units]
    return params

def updateParameters(self, params):
    if not params[1].altered and params[0].altered:
        spRef=arcpy.Describe(params[0].valueAsText).spatialReference
        if spRef.linearUnitName == "Foot_US":
            params[1].values = "0 Feet"
    return

I should also point out that this code defaults the distance to 0, where as the Buffer tool I'm imitating is able to set units without a distance.


texts: Generate a raster with a uniform value using QGis

texts:
I'm working on a project using the qgis-modeler and I want to generate a raster with a uniform value.

I'm using r.surf.random function of GRASS which needs min and max values but I want to have only one input.

So how can I generate it?


texts: Understanding the format: WKB from WKT and how to convert the first into the latter

texts:
I don't understand which kind of format I have my data. The column name is wkb_geom, so I supposed that data were in WKB format, but then I was checking around and I couldn't find example of it. Data are like this:

"0106000020E6100000010000000103000000010000007218000007000060B1D42B4010000060A372454007000060B1D42B40030000009D724540030000E0D5D42B40030000009D724540030000E0D5D42B40050000C08A7245400B000040FAD42B40050000C08A7245400B000040FAD42B40130000807B7245400B000040FAD4 (...)"

Is it in WKB or WKT format?? Second question, in case it's in WKB format, how can I convert it into WKT format? I was trying to follow this suggestion

How to convert WKB to WKT?

so the query is

UPDATE "ita_adm1"
SET wkb_geometry  = ST_GeomFromWKB("wkb_geometry",4326)

but it keeps saying that the function ST_GeomFromWKB doesn't exist.


texts: Inconsistent performance using the tiger geocoder.

texts:
I'm testing the geocoder results for one state, and I'm getting different results for directions that seem to be well formatted.

E.G.

Geocoding these addresses is taking around 100ms each:

280 FIRST AVENUE,MANHATTAN,NY,1000
214-17 PALMER DRIVE,QUEENS,NY,11697
130 LIBERTY STREET,MANHATTAN,NY,10006
345 PARK AVENUE,MANHATTAN,NY,10154

While for this group:

53 5TH AVENUE,MANHATTAN,NY,10033
810 7TH AVENUE,MANHATTAN,NY,10170
4602 13TH AVENUE,BROOKLYN,NY,1120

I don't find any clear pattern for to justify this difference, it seems to be happening for all the zip codes:


texts: How to change default password in GeoServer

texts:
I started to learn GeoServer to publish online maps.
I installed the software and now I would like to change the default password. At the moment is: geoserver
I tried to do it without result.
Can anyone explain me how to do this?


texts: How to read data from PostGIS in DotNet (VB.Net/C#)

texts:
I have stored my data in PostGIS (PostgrSQL). Now I want to read
the data in DotNet ArcEngine Application using VB.NET/C#.

Anyone know how to read the data from PostGIS in .Net?

Regards


texts: Drawing graph elements (nodes and edges) in Leaflet?

texts:
I need to build a simple online map editor and I'm intended to use leaflet (but can eventually use openlayers). How can I draw graph elements (nodes and edges) in leaflet (on top of a map). Specifically, I can draw linestrings (edges) and circles (nodes), but they are not correlated in any way, that is, the linestring is not an edge of the point since they are two different geometries. Hence, how can I draw graph elements, nodes and edges, where edges and nodes are associated.

Example:

Here, the circles should represent a node, while the polylines are the edges. However, the two geometries are not associated in any way, that is, I can't associate the edges with the nodes.

How do I draw a graph on top of such mach where I can add, retrieve, delete edges and nodes?


texts: inserting Cached .mxds into Flex viewer 2.5

texts:
I am trying to cache a number of operational layers in my flex viewer 2.5 application but after the layer type is changed from dynamic to tiled they are not visible in my viewer they are still visible in the Table of Contents widget. When I change them back to dynamic and reset they are still invisible. Is anyone familiar with this issue and have a solution?
Thanks,
jay


texts: carto api with javascript using online carto db. Argument of WHERE must be type boolean, not type geometry

texts:
I am using cartodb api to fetch the geojson from the cartodb data set using this code

var map = L.map('map',{ center: [-35.1393961, 149.8530818], zoom: 13});
var alpha=map.getBounds();

var lat0=alpha.getNorthEast().lat;

var long0=map.getBounds().getNorthEast().lng;

var lat1=map.getBounds().getSouthWest().lat;
var long1=map.getBounds().getSouthWest().lng;

// Add Tile Layer basemap
L.tileLayer('https://api.mapbox.com/styles/v1/mapbox/streets-v10/tiles/256/{z}/{x}/{y}?access_token=pk.eyJ1IjoiaGFzaGlyaWx5YXMiLCJhIjoiY2l4NzdzMmhjMDB0dDJ0bXN0NW5jc3cwayJ9.UNDyJQP92EC0ToFilX28lA', {
maxZoom: 18,
id: 'your.mapbox.project.id',
accessToken: 'your.mapbox.public.access.token'}).addTo(map);

// Global Variables
// Will go here
var coffeeShopLocations = null;

// Database Queries
// Get all coffee cafes from dataset
var sqlQuery = "SELECT * FROM zearch_nsw_19042016 WHERE the_geom && ST_MakeEnvelope("+long1+", "+lat1+", "+long0+", "+lat0+", 4326);";
function showAll(){
  $.getJSON("https://{username}.carto.com/api/v2/sql/?q="+sqlQuery+"&api_key={key}&format=geojson&callback=?", function(data) {

  console.log(data);

       });
};

// Run showAll function automatically when document loads
$( document ).ready(function() {
  showAll();
});

While executing the query on the cartodb it is returning correct result but here at first it was showing invalid geom object but when I outputed the geojson it is showing this error

  argument of WHERE must be type boolean, not type geometry

What is wrong?
The query works fine when implemented on the cartodb or postgis database. But when used with cartodb api it generates the error on geojson


texts:
I'm trying without success to find a solution for problem in ArcGIS and this is the reason that I'm sending this message.

I have one polygon shapefile with postal codes for South Africa (no Gap, no overlap). There are some small polygons that seem to be box codes inside big polygons (Street code) and I would like to identify or select them. The problem: I don’t have any attribute in the table to make a classification (Box code vs street code).

I'm looking for a tool or a method in ArcGIS to identify polygons inside polygons using just one shapefile.


texts: GDAL external overviews not recognized by ArcGIS

texts:
I am adding external overviews to a number of TIF files with GDAL as follows:

gdaladdo -ro file.tif 2 4 8 16 --config COMPRESS_OVERVIEW LZW

Then, when I add the tif file to ArcGis is still asking me if I want to build pyramids. It seems that the .ovr file is not fully recognized. Any suggestion?


texts: How to reproject WSF layer created from GeoJSON?

texts:
I have a WFS layer:

myVecLayer_disp= new OpenLayers.Layer.Vector("Редактируемый участок",{projection: new OpenLayers.Projection("EPSG:900913")});

And i have feature store:

        var protocol=new OpenLayers.Protocol.HTTP({
            url: url_servlet+"fullSearch.jsp",
            format: new OpenLayers.Format.GeoJSON({internalProjection: new OpenLayers.Projection("EPSG:4326"),externalProjection:new OpenLayers.Projection("EPSG:900913")})
    });
    searchStore = new GeoExt.data.FeatureStore({
        layer: myVecLayer_disp,
        fields: [
            {name: 'id', type: 'String'},
            {name: 'name_dispatcher', type: 'String'},
            {name: 'voltage', type: 'String'},
        ],
        proxy: new GeoExt.data.ProtocolProxy({
            protocol: protocol,
        }),
        autoLoad: true,
        layerOptions: function() {
            return {
                strategies: [new OpenLayers.Strategy.BBOX()],
                //projection:"EPSG:4326"
            };
        }
    });

As you see i use GeoJSON. In this geojson i have features in wgs84. And reproject it in mercator in protocol. Now i show this layer:

    app.mapPanel.map.addLayers(myVecLayer_disp);
myVecLayer_disp.visibility=true;

And see features near Africa but expect to see it within Russia. So Nothig reproject i think. But when i see layes's features:

    LINESTRING(0.00055807640304122 0.00050642147755534,0.00055807616948463
0.00050641930362077,0.00055807588202664 0.00050641664462617,0.00055807423809849
0.00050641332981391,0.00055807306129967 0.00050640979946183,0.00055807191147579
0.00050640644872202,0.00055807052802261 0.00050640231648177,0.00055806910873212
0.00050639813934481,0.00055806841697321 0.00050639609116937,0.00055806763547558
0.00050639388128181)

Ant this is same feature in geojson:

{"type":"LineString","coordinates":[[62.124781001642077,56.374581001212391],
[62.124755002241727,56.374338999945053],[62.124723002564679,56.374043002482722],
[62.124540001320156,56.373673997157425],[62.12440900067503,56.373281000617851],
[62.124281002866155,56.37290799870533],[62.124126997563053,56.372447999975407],
[62.123969002867341,56.371983002813522],[62.123891996618511,56.371755002191776],
[62.123805000700372,56.371508997142655]]}

What i doing wrong?


texts: How to merge waypoints or tracks into a single block?

texts:
When downloading routes from OverpassTurbo as KML files, I end up with two blocks: One is made of "<Placemark><LineString><coordinates>", and the other of "<Placemark><Point><coordinates>".

To merge those as a single route, I checked Viking and JOSM but didn't find how do it. I also tried GPSVisualizer, but this is what I get:

Does someone know of a way to merge ways or tracks into a single route?


texts: Extracting the grid coordinates in a raster map

texts:
I have a gridded raster map produced from the NASA precipitation data over a country. The pixel type of the map is in floating point with 32 Bit pixel depth. Therefore, it does not have any attribute table. Now I have two questions:

If I use the Int command in the ArcToolbox to convert the map to Integer point pixel type I lose lots of information of the gridded data. What is the solution to convert it to Integer point pixel type without losing so much information. The data are rainfall amount in 10km *10km grid points and the values are mainly different in regards to their decimal and removing that convert the map into just 4 categories!
In the original map with the floating point with 32 Bit pixel depth, how can I extract the 10km *10km coordinates? It is too many coordinates (above 15000 grids) and I need a solution to extract it. Converting to ASCII did not help, by the way.


texts:
I have GeoServer running locally, serving up a large table of vector points stored in Postgres. I have enabled a time dimension on GeoServer, from a column of type TIMESTAMP WITH TIMEZONE. from QGIS, I connect to the WFS. This all works great.

However I cannot use the feature subset of the layer properties in QGIS. The "Query Builder" button is greyed-out and the text box under "Feature subset" is not editable.

How can I issue a feature subset query for a WFS in QGIS? For instance, I want to only display points with a time property within a particular date interval. I can add filters to the query itself, but that is not easy to use compared to the query builder. Is this a GeoServer configuration issue? (I just started up a quick demo using Docker, and changed no GeoServer settings outside of those settings I had to choose when establishing my data store, etc.) I need to be able to subset my features, because the full dataset is around 500,000 records.

EDIT: I can open the attribute table and select features with an expression there, and these points are selected on the map. Its just the properties dialog does not allow me to subset.

EDIT: GeoServer 2.9 (snapshot), PostgreSQL 9.3 and QGIS 2.10 if these matter.

EDIT: Can confirm that what I want to do is possible only with the QGIS WFS version 2.0 plugin, but not the default add WFS dialog...


texts:
Please forgive my ignorance if I am asking a totally irrelevant question, I am still a GIS newbie.
I looked in Google and found a lot of links as of how to convert point clouds to DEMs.

I am trying to achieve exactly an opposite:
to create a point cloud from a DEM of the specific location and then
to create a 3D model of this location from this point cloud.

I didn't find one single link about it.

So is all my thinking wrong? Is there other way to achieve what I am trying to achieve?


texts: Import Pandas or any Python module into QGIS internal Python

texts:
(Also asked in stack overflow -
https://stackoverflow.com/questions/43244062/import-pandas-or-any-python-module-into-qgis-internal-python  )

I am developing a plugin for the current LTR version of QGIS (v2.14) using the Python installation included in it (Python v2.7).

This Python installation includes Numpy, but not Pandas. I am working with large datasets and performing large, slow operations on them, and Pandas would do some the work way more efficiently (faster) than using only Numpy - as I can judge when I run my scripts in a separate installation of Python 3 on my machine (on Windows 7). Unfortunately Pandas is not included in the Python 2.7 installation of QGIS.

Installing Pandas manually is an option if the plugin was to be used by 1-5 users, but I plan to distribute to a large group of users (using Win7 and Win10) so using Pandas would require an automatic way to include Pandas in the plugin somehow, or automatically install it for the internal Python QGIS installation.

Is there any automatic way (a Python script) that I could call at the beginning of my script to:

1) just call Pandas from a folder and import it without errors or

2a) check if Pandas is installed in QGIS-Python

2b) if not, install it in the internal QGIS Python version (without messing too much with environment variables for the PCs, as some of them have alternate versions of Python in use).

So something like

# My Python plugin script starts here

import numpy as np
import math
import ...

functionPandasPresent():
    #function to check if Pandas is present here
    #result = True or False, so for example
    result = False
    return result

pandas_present = functionPandasPresent()

if pandas_present is False:
    # function to install Pandas here
    # looking for whatever dir QGIS-Python was installed in
    # installing it just for that Python installation

import pandas as pd

# The rest of my script continues here...

With a script like this, not only Pandas but any other additional package could be used in custom QGIS plugins - so useful for any plugin developer I guess.


texts: How to remove a QgsMapCanvasLayer from QgsMapCanvas with PyQGIS || C++?

texts:
I try to find a way to remove a QgsMapCanvasLayer from the QgsMapCanvas. At the moment I only can delete the MapLayer from the QgsMapLayerRegistry with:
    QgsMapRegistry::instance()->removeMapLayer(LayerID);

The MapLayer vanishes from the QgsMapCanvas, but: macanvas->layerCount();
results in the same amount as before. My Question is:
Do I have to create a new/updated QList<QgsMapCanvasLayer> as layerSet in order to the updated QgsMapLayerRegistry, or is there a method to remove the deleted QgsMapLayer directly from the QgsMapCanvas
Or maybe an easier question. Can I cast QList<QgsMapLayer*> to QList<QgsMapCanvasLayer>? Or do I understand the concept wrong? Kindly, olleonkel


texts: Find minimum number of polygons connecting two points

texts:
Given GIS data of polygons and two points, is there a way to find the minimum number of polygons you need to traverse to walk between the two points? Even an approximate approach would be fine for my purposes if the problem cannot always be solved exactly.

The actual application I am using this for uses tax parcels and pipelines, with the goal of finding minimum number of parcels I need to traverse to lay a pipeline between two points -- so thousands of polygons and thousands of pairs of points. But to provide an example here is some code in R using the sf package. I load North Carolina counties and put two points on it. The straight line path between the two points goes through 5 counties. But if you look at the map you can see there's a path that goes through only 4 counties. I want an approach to calculate that 4 in this case, and more generally.

# install.packages(c('sf','dplyr')) # if needed
library(sf)
library(dplyr)

# use NC counties polygons
polys <- read_sf(system.file("shape/nc.shp", package = "sf"))

# make two points, goal is to find the minimum # of counties to get between them
points <- data.frame(name=c('start','end'),LAT=c(34.623858,35.257017),LONG=c(-78.851954,-78.237195))
points <- st_as_sf(points, coords = c("LONG", "LAT"), crs = 4326)

plot(polys$geometry)
plot(points$geometry, add=TRUE,col='blue')
plot(st_cast(summarize(points,name=first(name)),'LINESTRING'),add=TRUE,col='blue')


texts: Checking if geometry in one table exists in another using PostGIS?

texts:
I am trying to find the number of animals each state has. I first did this query:

    COUNT(*)  from animals c, locations a where st_contains(a.geom,c.the_geom);

But for some reason I get a syntax error. Can someone guide me in the right direction?

Below are examples of my two tables.

Animals table

      animal_type     city area     the_geom
       9            243582        82F4100000000F0BB0D41

       4           11234         6C892E410000000040170

       3           243582        82F4100000000F0BB0D41

      12          11234         6C892E410000000040170

Location table

     State        gid                  geom

     WV           3                    82F4100000000F0BB0D41

     TX             2                    6C892E410000000040170


texts: Ingesting LARGE CSVs to file geodatabase to create feature class and then clipping?

texts:
I am trying to convert 30 .csv files that have ~50,000,000 records each, and are just over 3gb each, to point feature classes within a FGDB.

If I do this one by one by hand, I generally right click in ArcCatalog > Create Feature Class > From XY Table--this takes a while but works.  Unfortunately I cannot call that tool from arcpy when trying to automate this--doesn't exist to my knowledge.  I wrote the script below to automate the ingest.  It works on small files, but fails on these large files.

Any advice?

Willing to explore other tools.

import csv
from itertools import islice
import os
import arcpy

arcpy.env.overwriteOutput = True
fgdb = r'C:\_L\Cols.gdb'
clip_features = r'C:\_L\MS.gdb\Polys'

for subdir, dirs, files in os.walk(r'C:\_L\Samples'):
    for file in files:
        if file.endswith(".csv"):
            fullFilePath = os.path.join(subdir, file)
            print fullFilePath
            fileName = os.path.basename(fullFilePath) #filename with extension
            fileNameNoExt = os.path.splitext(fileName)[0] #filename without extension

            arcpy.MakeXYEventLayer_management(fullFilePath, 'Field3', 'Field2', fileNameNoExt)
            arcpy.CopyFeatures_management(fileNameNoExt, fgdb + "\\" + fileNameNoExt)
            #arcpy.Clip_analysis(fgdb + "\\" + fileNameNoExt, clip_features, fgdb +"\\"+ fileNameNoExt + "_clip")
            print fullFilePath


texts: Applying Winner Takes All method to raster data in ArcGIS Desktop?

texts:
I have a classified raster dataset (A). I've transformed it to polygon type data and then have applied a fishnet to it.(B) Now i want to assign a single class to a single grid (polygon type i need) with the "winner takes all" method (C).

Am i doing right with changing data type from raster to polygon?


texts: Using ogrinfo to list all projection source of a database

texts:
I have to reproject all the layers from a database, but some of them are undefined and some are defined in a wrong projection. I would like to sort all those layers in a text file with their source projection and their path. I tried:

For /R %f in (*.tab) do ogrinfo -al -so %~dpnxf > list.txt

The fact is only the last layer return a result in the list.txt. Do you know why? Is there a more appropriate method to overcome that projection problem?


texts: Zonal histogram from 3 band RGB soil raster

texts:
I have soil map with 3 RGB bands and I need a zonal histogram for polygons from this raster (how many pixels of each soil the polygon contains).

First I tried to sieve the raster and then get zonal histogram, but I found that sieve tool by default uses first band (red in my case). So different colors with the same value in the red band get the same value after the sieve. I tried this with another band but result was the same. Also I tried to create raster layer with raster calculator (just sum over values from all 3 bands).

Can anybody recommend any solution using QGIS?


texts: Importing DBF files into Microsoft Excel

texts:
I'm working through the book "GIS Tutorial 1" and it wants me to open up a DBF file in Microsoft Excel. But when I try to open it, I'm running into some issues. First thing that seems weird is although the extension is .dbf, it's listed as an XML document.

When I try to open it up in Excel, instead of getting the neat table I expected I got a bunch of additional information.

I opened up the file as listed as 'OpenOffice.org ...' and it mostly gave me what I wanted, but some of the text is left as ####.

How can I get a neat table like in picture 3, but using a DBF file? Also how can I get it to include all data, and not block some out with #####?


texts: How to serve a large vector dataset on GeoServer

texts:
I have a GeoServer on a VPS (6 cores, 24GB RAM).
I created PostGIS store with a dataset of around 12GB (4.4 M polygons).
From store I imported 10 layers (each represents a region of a country). It takes 2-3 seconds to load a WMS tile on average from these layers.
I want to serve all of these layers as a single service, but when I import all of them into a single Layer Group, it takes more than 60 seconds to render a single WMS tile, so I receive a TIMEOUT error.

Example:
Loading a WMS tile from layer A (bbox 0-10) takes 2 seconds.
Loading a WMS tile from layer B (bbox 10-20) takes 2 seconds.
...
Loading a WMS tile from layer J (bbox 90-100) takes 2 seconds.
Loading a WMS tile from Layer Group A-J (bbox 0-100) takes 60+ seconds. (Even if this tile is within layer A bbox only).

I expected this Layer Group to load a tile in 4s~ max. There should be some sort of major performance loss that I don't know about. What is the bottleneck here? Maybe I should use different technology for this?


texts: Mapbox GL control visibility on multiple labels using zoom attribute

texts:
I have a point dataset that is intended only as labels. Each label has a pre-calculated minZoom attribute to one decimal place, e.g. 6.1. While there are labels that share the same value, they vary considerably. I've figured out how to programmatically add a single label. It's workable, I suppose, to add a layer for each label. But performance may suffer. Is there a way, either via Mapbox Studio or programmatically (for the iOS SDK) to set the camera stops individually for each label but keep them within a single layer?

func mapView(_ mapView: MGLMapView, didFinishLoading style: MGLStyle)
{
    // Try adding point with label
    let coordinates: CLLocationCoordinate2D = CLLocationCoordinate2D(latitude: 49, longitude: -114)

    let testPt = MGLPointFeature()
    testPt.coordinate = coordinates
    testPt.title = "Testing label"
    testPt.attributes = ["title": testPt.title as Any]

    let src = MGLShapeSource(identifier: "testPtID", features: [testPt], options: nil)

    mapView.style?.addSource(src)

    let symbols = MGLSymbolStyleLayer(identifier: "testPtSym", source: src)

    let color = UIColor(red: 0.08, green: 0.75, blue: 0.96, alpha: 1.0)

    symbols.textColor = MGLStyleValue(rawValue: color)
    symbols.text = MGLStyleValue(rawValue: "{title}")
    symbols.textFontSize = MGLStyleValue(rawValue: 24)
    symbols.textOpacity = MGLStyleValue(interpolationMode: .exponential,
                                        cameraStops: [
                                            6: MGLStyleValue(rawValue: 0),
                                            6.1: MGLStyleValue(rawValue: 1)
        ], options: nil
    )

    mapView.style?.addLayer(symbols)
}


texts: Setting up Geoprocessing Model with ArcGIS for Utilities CCTV Processor?

texts:
I have tried to find any already available data on this issue but have come up empty handed.

I am using ArcMap 10.3.1 and MS SQL Server in a versioned environment, for inspections we are using a third party CCTV inspection software called pipetech.

I am attempting to check out CCTV Processor for ArcGIS Utilities (http://solutions.arcgis.com/utilities/water/help/cctv-processor/) I have downloaded all the data and models from the ESRI page and I have begun to change the parameters to my MS Access CCTV dump, etc...  The cctv software (3rd party) we use for inspections writes its data to a MS Access database. The model from ESRI will access data from this DB and use it to tag pipes and create a NASSCO rating for the section (at least this is what I am assuming)

The model that they prebuilt for the tool is very complicated to say the least and I am having a problem figuring it out, when I change certain things I seem to completely crash the model.

Does anyone have a tutorial or help file that would help explain how the model works and which components of the model need to be changed?


texts:
Simply I want to upload a tiny csv file to a buffer Geoprocessing service published through ArcGIS Server 10.4.1. The problem ,however, is that I can't access any file in my server or an external server no matter what method that I use. I have tried uploading the file directly to the GP tool -I can upload it from the service directory- , using request function and even defining csvURL that is actually working in an example published by ESRI, sometimes even those examples when I download their sources everything works fine except the layers csv,kml,... I can't call them.
I don't know what is blocking my requests is it the firewall or the proxy?
I am still new to this field.


texts: Live connection between ArGIS Online and Google Docs Spreadsheet?

texts:
I want to make a join in ArcGIS Online: From a spreadsheet with names (visitors to an exhibition) to polygons (containing that name).
I succeed to make that in ArcGIS Online: but there is no live connection.

Has anyone already done this?

Do I really need Zapier to make that connection?

We have an organisational account on ArcGIS Online

I have already tried several things:
add an item from cloud drive
within an item: add a layer from the web (.csv spreadsheet)

There is no need to change that spreadsheet from within ArcGIS Online.


texts: Saving Python objects into memory for repeated use of Python toolbox?

texts:
I have made a Python toolbox which processes relatively large polyline feature class (millions of features), builds a network graph from it as dictionary of lists, and depending on user selected features in map finds shortest paths between selected features and creates a new selection of these shortest paths. This tool is used repeatedly in one of our workflows. I would love to keep the "network graph" as dictionary of lists in memory. I am aware of in_memory workspace and possibility to convert my dictionary into table and store it like that there.

I wonder if there is some way to save directly Python objects into generally accessible space (variable) in memory in ArcGIS Pro from Python toolbox?

Something like in_memory workspace but directly for Python objects (variables).
In my particular case so I can check if this dictionary of lists exists and potentially reuse it if object meets the user settings and other requirements, when the tool is rerun?


texts:
I want something similar to OpenLayers.Format.OSM for leaflet. i know how to use a geojson file and i have found support for GPX and WFS in plugins, but not for OSM.

Simply saying i want to make this map using leaflet and overpass.


texts: Using multipatch and surface in ArcGIS 3D Analyst?

texts:
I want to get 3D building. I have DEM surface and 3D feature of building roof. I created multipatch class by extruded building roof with '-' value.

Is it possible to create building from this multipatch and DEM surface ?

Interpolate Polygon To Multipatch (3D Analyst) create only footprnt of building's roof on TIN surface. Output is Multipatch but it is only plane. How can I create Multipatch between Roof's polygon and footprint's Multipatch ?


texts: Update Service config file with ArcIms java connector

texts:
We are using arcims java connector to create, update, delete services. it works fine for create and delete. However when we try to update the config file of service (ie : mxd file) the new mxd file is not replaced in the ArcImsSite.sez file. Thus when we restart ArcIms, the service is created with the wrong mxd file.

here is a snippet of what we do to update the service :

ConnectionProxy connectorProxy = new ConnectionProxy();
connectorProxy.setConnectionType(ConnectionProxy.TCP_ADMIN);
connectorProxy.setHost(arcimsHost);
connectorProxy.setPort(arcimsPort);
connectorProxy.setDisplayMessages(false);

configfile = new FileInputStream(ConstantManager.ENV_DIRECTORY+"configuration.properties");
properties.load(configfile);

int ping = connectorProxy.ping(false);
Site site = new Site();

if (ping == ConnectionProxy.PING_FIRST_LOGIN) {
    site.setSiteUser(properties.getProperty("ADMIN_AIMS"),
    properties.getProperty("PASSWORD_AIMS"),
    connectorProxy);
} else if (ping == ConnectionProxy.PING_OK) {
    connectorProxy.setUsername(properties.getProperty("ADMIN_AIMS"));
    connectorProxy.setPassword(properties.getProperty("PASSWORD_AIMS"));                            site.setSiteUser(properties.getProperty("USER_AIMS"), properties.getProperty("PASSWORD_AIMS"),                          connectorProxy);
} else if (ping == ConnectionProxy.PING_FAIL) {
LoggerManager.logError(LOG, ConstantManager.BUSINESS_TYPE, null);
throw new BusinessException(AdminErrorMessage.ARCIMS_NOT_AVAILABLE, null);
}

ServiceCollection serviceCollection = ServiceCollection.getServices(connectorProxy);

for (int i = 0; i < serviceCollection.size(); i++) {
  Service aservice = serviceCollection.getService(i);
  // delete the old service
  if (aservice.getName().equals(mapBefore.getMapName())) {
   LOG.debug("\t-->INPUT BUSINESS  LAYER on  updateMap METHOD : suppression du service ");
   serviceCollection.removeService(aservice);

aservice.removeService(connectorProxy);
  }
}

//Here i try to save the site config with site.save(connectorProxy). However the service is //deleted from ArcImsSite.sez file by doing this

//Recreate the connexion
ConnectionProxy connectorProxy = new ConnectionProxy();
connectorProxy.setConnectionType(ConnectionProxy.TCP_ADMIN);
connectorProxy.setHost(arcimsHost);
connectorProxy.setPort(arcimsPort);
connectorProxy.setDisplayMessages(false);
int ping = connectorProxy.ping(false);
Site site = new Site();

if (ping == ConnectionProxy.PING_FIRST_LOGIN) {
    site.setSiteUser(properties.getProperty(USER_AIMS"),
    properties.getProperty("PASSWORD_AIMS"),
    connectorProxy);
} else if (ping == ConnectionProxy.PING_OK) {
    connectorProxy.setUsername(properties.getProperty("USER_AIMS"));
    connectorProxy.setPassword(properties.getProperty("PASSWORD_AIMS"));                            site.setSiteUser(properties.getProperty("USER_AIMS"), properties.getProperty("PASSWORD_AIMS"),                          connectorProxy);
} else if (ping == ConnectionProxy.PING_FAIL) {
LoggerManager.logError(LOG, ConstantManager.BUSINESS_TYPE, null);
throw new BusinessException(AdminErrorMessage.ARCIMS_NOT_AVAILABLE, null);
}

//Create the new service. Actually same service as previous but with new mxd file.
Service myservice = new Service();
//.... adding required information to new service
myservice.setConfigFile(mxdFilePath);
myservice.setConfigFileLength(buff.length);
myservice.setConfigContents(bais, true);
//....

//Save current service config
if (!site.save(connectorProxy)) {
throw new BusinessException(AdminErrorMessage.SERVICE_CONFIG_CANT_BE_SAVED, null);
}

Does anyone have an idea of what is wrong ?


texts: GIS file format for space and time data?

texts:
I need some advice with GIS file format choice.
Here is my problem:
I have about 100 stations.
At each station water temperature is measured daily for 1 year
I want to load this data to QGIS and be able to quickly display a water temperature map for a user-selected day.
I came up with following file structure: A shapefile with lat, lon, time, value columns. For each combination of lat, lon, time there is a new point shape. So for each station I have 365 "overlapping shapes". Then I can use the "query" in QGIS to filter out shapes that have the time attribute equal to the selected day.
This is OK for 365 days but becomes little inefficient for several years (shapefile becomes huge and QGIS displays it very slowly)

Is there any other GIS format (like Spatialite or Postgis) that I could use to store the time-space data more efficiently and still be able to visualize it in QGIS without extra programming?


texts: Plotting heatmap on QGIS using external CSV values rather than density of lat/long points

texts:
I'm trying to plot a heatmap in QGIS.
The lat/long coordinates I have are arranged into 4 parallel lines (4 survey transects), with a third column in the CSV file comprising a value ranging from 1 to 49 (a count of the number of species per point).
I have managed to produce a map that simply plots these coordinate points on a basemap, and then apply a heatmap based on density of the points. These coordinates are however taken at regular intervals, and therefore the heatmap is pretty homogenous.
How therefore would I alter it to draw the heatmap based upon the third column of values?


texts: Convert tile naming from crs3912 to crs3857

texts:
I adapted this solution that projects crs3912 map to leaflet(crs 3857 if I'm correct) to node.js app that allows me to load this map into applications like Global Mapper, OruxMaps or such. However, zoom and map placement is off place.

I found that the tile https://c.tile.openstreetmap.org/15/17795/11679.png geographically corresponds to this geopedia tile, which when converted to decimal would be

  8/167/162 (z/x/y).

This map goes only till the zoom level 9 where the jpg files are already pixelated, so last usable level is 8.

My question is, simplified, when z/x/y coordinates requested from navigation app are for example 15/17795/11679 how do i get this tile, how to adapt coordinate system and zoom for this to usable in mentioned applications.


texts:
Below is my code which i am using for modifying features and drawing new features.

saveStrategy = new OpenLayers.Strategy.Save();
saveStrategy.events.register('success', null, ChangesSuccess);
saveStrategy.events.register('fail', null, ChangesFailed);

wfs = new OpenLayers.Layer.Vector("Edit Vagmitt Polyline", {

    styleMap: StyleMap,
    eventListeners: {
        "featuresadded": ZoomToData
    },
    strategies: [new OpenLayers.Strategy.Fixed(), saveStrategy],
    protocol: new OpenLayers.Protocol.WFS({
        url: GEOSERVER_WFS,
        featureType: "new_Vagmitt_polyline",
        featureNS: "http://www.openplans.org/workspace",
        geometryName: "the_geom",
        srsName: "EPSG:900913",
        version: "1.1.0",
        extractAttributes: true,
        isBaseLayer: false,
        visibility: true
    })

});

new OpenLayers.Control.DrawFeature(wfs, OpenLayers.Handler.Path, {displayClass: 'lineButton', title: 'Draw new line feature', handlerOptions: {style: sty}});

new OpenLayers.Control.ModifyFeature(wfs, {title: 'Edit line feature'});
The code works fine for modifyFeature control, but its gives error on save strategy when i draw a new feature.

The changes are also saved back to Geoserver when i modify the feature (triggers ChangesSuccess). but in case of DrawFeature control, saveStrategy triggers ChangesFailed function.

I could not understand, why ModifyFeature works but DrawFeature does't work.
Can anybody explain, where is the problem ?


texts: How to create an AAT table in ArcGIS Desktop?

texts:
A long time ago I used to create aat.dbf tables in ArcInfo by using the simple command "BUILD."
I'm staring at ArcGIS 9.3 and I have no clue about how to perform this same, simple task!

Is there a quick way to do this?


texts: Logging geoserver's request logs into postgres

texts:
I'm trying to force geoserver to write request.log(s) into a postgis database. Unfortunately, geoserver's documentation remains silent about logging into database at all.

Searching the web I found one answer on mailing list for similar question. I didn't try this, but that's because I've got feeling that's not thing I'm looking for, because this solution seems to allow to save to database logs of geoserver's activity not incoming request what is my point...

I also found geoserver's plugin called monitor, which could be thing I'm looking for, (again) unfortunately in v 2.3 adding it brakes answers for Request=GetFeatureInfo, so it's useless for me.

Finally, there's always a solution to use PostgreSQL COPY to import data from logfiles into database, but I think there have to be more elegant way.

The point is to log requests into database, so after little post-processing on request strings, I'll be able to write simple reports of IPs using my server or most commonly downloaded areas.

More specified question is: where can I find settings for request logs and how to reconfigure it to put data into database?

Can anybody help?


texts:
I have several annotations in my map that disturb.

Is it possible to hide them somehow and show them later when I need them?


texts: How to convert UTM coordinates into Decimal Degree coordinates?

texts:
I work on QGis 2.0.1 - Dufour on Macintosh OS X Mavericks.

I have UTM coordinates (XXXX and YYYYY) in a table (with one column of X and an other one of Y). I need this coordinates in Decimal Degree. Could you tell me the detailed protocol to follow please ?


texts: Seeking Vehicle Routing Problem and Location-Allocation Tool in QGIS?

texts:
I was wondering if there was a tool in QGIS that is similar to ArcGIS Network Analyst's Vehicle Routing Problem and Location-Allocation tools?

I could not find any info about this in the internet.


texts: Deleted Topology and Feature Classes are still locked within Geodatabase

texts: I deleted my topology within a geodatabase rather than disconnecting the feature classes inside of it. Doing this appeared to have deleted all 4 feature classes that were within the topology entirely from the geodatabase. I was able to find a backup of the geodatabase, but now when I try to export the replacement feature class into the original geodatabase, I am getting the error that a feature class with that same name already exists. Except it does not appear within the geodatabase and I have checked in both ArcCatalog and ArcMap. How do I locate these feature classes that are deleted but still somehow locked within the geodatabase?


texts: Storing GDAL VRT sources inside QGIS project files (.qgs, .qgz, .gpkg)?

texts:
After I managed to store SVG symbols in Geopackage BLOBs (Read (and write) SVG from Geopackage with QGIS) another idea was growing in my mind: Maybe we can store simple GDAL VRT sources in GPKG BLOBs or *.qgz files, to bind these datasources to a QGIS project?
Here's a minimal VRT which shows the city of Vienna clipped from a public WMTS (https://basemap.at):
<VRTDataset rasterXSize="293885" rasterYSize="227105" subClass="VRTWarpedDataset">
  <SRS dataAxisToSRSAxisMapping="1,2">EPSG:3857</SRS>
  <GeoTransform>  1.8013609905169420e+06,  1.4929107087160037e-01,  0.0000000000000000e+00,  6.1606686197213698e+06,  0.0000000000000000e+00, -1.4929107087160004e-01</GeoTransform>
  <VRTRasterBand dataType="Byte" band="1" subClass="VRTWarpedRasterBand">
    <ColorInterp>Red</ColorInterp>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="2" subClass="VRTWarpedRasterBand">
    <ColorInterp>Green</ColorInterp>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="3" subClass="VRTWarpedRasterBand">
    <ColorInterp>Blue</ColorInterp>
  </VRTRasterBand>
  <VRTRasterBand dataType="Byte" band="4" subClass="VRTWarpedRasterBand">
    <ColorInterp>Alpha</ColorInterp>
  </VRTRasterBand>
  <BlockXSize>512</BlockXSize>
  <BlockYSize>128</BlockYSize>
  <GDALWarpOptions>
    <SourceDataset relativeToVRT="0">WMTS:https://www.basemap.at/wmts/1.0.0/WMTSCapabilities.xml,layer=bmaporthofoto30cm</SourceDataset>
    <Transformer>
      <ApproxTransformer>
        <MaxError>0.125</MaxError>
        <BaseTransformer>
          <GenImgProjTransformer>
            <SrcGeoTransform>977649.895640023053,0.149291070871599985,0,6281289.98263385892,0,-0.149291070871599985</SrcGeoTransform>
            <SrcInvGeoTransform>-6548616.00182950962,6.69832424780491387,0,42074116.9981704801,0,-6.69832424780491387</SrcInvGeoTransform>
            <DstGeoTransform>1801360.99051694199,0.149291070871600373,0,6160668.61972136982,0,-0.14929107087160004</DstGeoTransform>
            <DstInvGeoTransform>-12066100.0018294789,6.6983242478048961,0,41266155.9981704652,0,-6.6983242478049112</DstInvGeoTransform>
          </GenImgProjTransformer>
        </BaseTransformer>
      </ApproxTransformer>
    </Transformer>
    <BandList>
      <BandMapping src="1" dst="1" />
      <BandMapping src="2" dst="2" />
      <BandMapping src="3" dst="3" />
    </BandList>
    <SrcAlphaBand>4</SrcAlphaBand>
    <DstAlphaBand>4</DstAlphaBand>
    <Cutline>MULTIPOLYGON (((5501234.2650394 1052508.13225313,5501234.2650394 791668.70254492,5834964.1589386 791668.70254492,5834964.1589386 1052508.13225313,5501234.2650394 1052508.13225313)))</Cutline>
  </GDALWarpOptions>
</VRTDataset>

We could store the VRT file in a QGIS *.QGZ project file and retrieve it via /VSIZIP/, but unfortunately it will be deleted after a project save. (Maybe the core devs could define a QGZ folder where we can store some custom data which won't get lost?)
Geopackage BLOBs are no option either, because I couldn't find a way to open a data stream (i.e. base64 encoded) as VRT data source. It's possible to retrieve data from BLOBs using Python and open it in QGIS from a temp file, but I would prefer a more direct solution.
GDAL has so many Virtual File System handlers, maybe we could use VSIMEM or VSISTDIN to read the data from a Geopackage BLOB.
Does anyone have a brilliant idea to solve this tricky task without using too much Python code?
(@Even Rouault ... something for you?)


texts:
Using OpenLayers3 I have observed an incorrect alignment of the text when styling the features of the Cluster source. Also the baseline of the text appears to be ignored.

The below code sample (or the jsfiddle) presents how I am styling the features of the Cluster source. In short my goal is to present a single feature differently to a clustering of features.

    var aisLayer = new ol.layer.Vector({
        name: 'AIS Layer',
        source: clusterSource,
        style: function (feature, resolution) {
            var size = feature.get('features').length;
            var style;
            if (size === 1) {
                var theFeature = feature.get('features')[0];
                style = [
                new ol.style.Style({
                    image: new ol.style.Circle({
                        radius: 15,
                        stroke: new ol.style.Stroke({
                            color: '#000000'
                        }),
                        fill: new ol.style.Fill({
                            color: '#ffffff'
                        })
                    }),
                    text: new ol.style.Text({
                        textAlign: "Start",
                        textBaseline: "Middle",
                        font: 'Normal 12px Arial',
                        text: theFeature.get('name'),
                        fill: new ol.style.Fill({
                            color: '#ffa500'
                        }),
                        stroke: new ol.style.Stroke({
                            color: '#000000',
                            width: 3
                        }),
                        offsetX: 0,
                        offsetY: 0,
                        rotation: 0
                    })
                })];
            } else {
                style = styleCache[size];
                if (!style) {
                    style = [
                    new ol.style.Style({
                        image: new ol.style.Circle({
                            radius: 12,
                            stroke: new ol.style.Stroke({
                                color: '#ffcc33'
                            }),
                            fill: new ol.style.Fill({
                                color: '#000000'
                            })
                        }),
                        text: new ol.style.Text({
                            textAlign: "Center",
                            textBaseline: "Middle",
                            font: 'Normal 12px Arial',
                            text: size.toString(),
                            fill: new ol.style.Fill({
                                color: '#ffcc33'
                            }),
                            stroke: new ol.style.Stroke({
                                color: '#000000',
                                width: 1
                            }),
                            offsetX: 0,
                            offsetY: 0,
                            rotation: 0
                        })
                    })];
                    styleCache[size] = style;
                }
            }
            return style;
        }
    });

The problem: The text is not centred for the clustered features (black circles), as if it still has the alignment setting used by the not clustered features (white circles). Also the baseline of the text appears to be ignored. To see the problem in action go to the jsfiddle.


texts: Modify a Python script to one that uses RasterDomain_3d fuctionality - and writes resulting shapefiles it to one

texts:
As I am not very experienced with Python I need your help. I managed to modify a Script (that has been posted here Get boundary of raster image as polygon in ArcGIS Desktop? by https://gis.stackexchange.com/users/2043/jeb)

import arcpy,os

InFolder = arcpy.GetParameterAsText(0)
Dest=arcpy.GetParameterAsText(1)

arcpy.env.workspace=InFolder
#The raster datasets in the input workspace
in_raster_datasets = arcpy.ListRasters()

arcpy.CreateFeatureclass_management(os.path.dirname(Dest),os.path.basename(Dest),"POLYGON")
arcpy.AddField_management(Dest,"RasterName", "String","","",250)
arcpy.AddField_management(Dest,"RasterPath", "String","","",250)

cursor = arcpy.InsertCursor(Dest)
point = arcpy.Point()
array = arcpy.Array()
corners = ["lowerLeft", "lowerRight", "upperRight", "upperLeft"]
for Ras in in_raster_datasets:
    feat = cursor.newRow()
    r = arcpy.Raster(Ras)
    for corner in corners:
        point.X = getattr(r.extent, "%s" % corner).X
        point.Y = getattr(r.extent, "%s" % corner).Y
        array.add(point)
    array.add(array.getObject(0))
    polygon = arcpy.Polygon(array)
    feat.shape = polygon
    feat.setValue("RasterName", Ras)
    feat.setValue("RasterPath", InFolder + "\\" + Ras)
    cursor.insertRow(feat)
    array.removeAll()
del feat
del cursor

that writes the Extents of a List of rasters that are in one folder to a Shapefile and adds the name of the raster to the attribute table.
I would need the same function, but instead of the Extent of each Raster with the Raster Domain - to have the exact frame of each extracted georeferenced video still.

Because I do not completely understand the single steps... I do not know where to incorporate the arcpy.RasterDomain_3d...


texts: Copy polygons to another layer with different attribute

texts:
I am working in QGIS 3.0.2-Girona.

I am digitizing images in which there are rocks that can be encrusted, totally or partially, by algae. I want to have a polygon for the rock and one for the algae. In particular I can have 3 cases:

the polygon of a rock and the one of the algae coincide
the polygon of a rock could potentially contain more than one polygon of
the algae
the polygon of a rock could contain NO algae polygon(s) at all.

For this reason I thought that one effective way to go is to have 2 polygon vectors open in edit mode: 1 for rock and 1 for algae. When digitizing a rock with algae, I digitize the rock first and then copy paste the Geometry only (i.e. without attribute) to the algae vector layer. Repeat this process whenever I find a similar case while digitizing.

My problem is that selecting a feature and copy/paste it to a different layer doesn't seem to work: although the 2 layers are in edit mode, when I paste the geometry to the target layer I can see the pasted feature in its attribute table (with just the ID value, that's OK) but it doesn't appear on the map. If I select the feature in the target layer's attribute table, and click "zoom to selection", I am shown a white part of the map (I don't even know where) clearly not in the same position as the copied polygon.

Am I doing something wrong? Does anybody have a strategy to achieve this task?


texts: Making a polygon from list of xy coordinates using Shapely

texts:
I am new to Python (moving over from Matlab) and am trying check if a set of points are inside a polygon using Shapely as I am not able to install GeoPandas (tried for over a week now and am tired of fighting with Anaconda).
Background:
I have a folder of las files with coordinates in the filename and what I want to do is:

read the filename
parse it to x and y coordinates
check if this xy coordinate is inside of a polygon

I am able to do the first 2 steps but am failing on step 3. My problem is that I have a set of polygon coordinates (xy) in a .csv file that are being read into a list:
Sample of the data :
156922.3097 6559238.263
157424.9655 6559462.937
157837.7369 6559395.169
158078.2783 6559284.668
158449.978 6558972.523

The code I am using to read the csv file:
with open("xy_polygon.csv", 'r') as f_input:
    csv_input = csv.reader(f_input, delimiter=' ', skipinitialspace=True)
    xv = []
    yv = []
    for cols in csv_input:
        xv.append(float(cols[0]))
        yv.append(float(cols[1]))

and I am able to add the points by typing out something like this:
# x is a list of x-coordinates, y is list of y-coordinates
p1 = Polygon([[xv[0], yv[0]], [xv[1], yv[1]], [xv[2], yv[2]]])

But I would like to use both lists in the Polygon function, for example:
p1 = Polygon([[xv, yv]])

But I get the error:
ValueError: A LinearRing must have at least 3 coordinate tuples

I also tried:
p1 = Polygon([[xv[0], yv[0]], [xv[1], yv[1]], [xv[2], yv[2]]])
p1.append(Polygon([[xv[3], yv[3]], [xv[4], yv[4]], [xv[5], yv[5]]]))

which returned the following, so I am clearly doing this incorrectly.
'Polygon' object has no attribute 'append'

Can someone provide some tips as to how I could achieve this?
Edit: fixed error in original question/code


texts: heatmap plugin slow/unresponsive in QGIS 1.8

texts:
I am trying to use the new heatmap plugin to create a heatmap from about 10k points. When I enter the parameters and start the process , QGIS becomes unresponsive. I don't get a progress bar or a dialog. If I let it be, it generates the heatmap after an hour or so. On the same system, GRASS's v.kernel creates the heatmap in under a minute. I tried with a few different parameters and reducing number of points but it isn't much better. What am I doing wrong? Have others seen such behavior? I am using QGIS 1.8 on Windows 7.

The data I am trying to work with is at http://data.london.gov.uk/datastore/package/policeuk-crime-data


texts: Creating Directional Flow Arrows for sewer lines?

texts: I am working for the Town of Yadkinville, NC and would like to know how to create a layer of directional flow arrows for sewer lines to be able to turn them on and off at will?


texts: Accessibility analysis using isochrones and network analysis tools QGIS 3.16

texts:
So I made an isochrone using ORS Tools plugin, as well as using the Network analysis toolbox to analyze the service area of the bus stop. Ultimately I want to provide a statement of X% Area covered from the bus stops. So to do that I need the area of the service area of the bus stops divided by the area of the district (The yellow shapefile layer). It'll be easier to calculate the area if I use isochrones from the ORS Tools plugin, but it went beyond the district itself making the Service Area/District Area formula questionable, meanwhile, the network analysis toolbox is perfectly bounded inside the district layer but it's impossible to calculate its area since it's counted as a line vector.
So is there any method to keep the isochrone under the boundary of the district shapefile layer? Or is there any method to turn the service area lines into polygons to calculate its area?


texts: How to interpolate a lat/long point between two others at a given distance

texts:
I'm trying to solve the problem of interpolating a point at a given distance (in meters) between two other lat/long points. (e.g. P1 (50,10), P2(50.01, 9.990), I want to interploate P3 at a distance "X" on the shortest route between P1 and P2, assuming that X is always smaller than the distance between P1 and P2)

Whilst at first this seemed kind of easy I couldn't for the life of me find an easy to implementation for this (e.g. a few lines of code and a few external constants, as to be easy to debug and quick to run).

Currently I'm using the formula:

val x : Double = RadiusOfEarth*Math.toRadians(lon)*CentralLatCos
val y : Double = RadiusOfEarth*Math.toRadians(lat)

to project to 2d cartesian, I interpolate the point in cartesian coordinates (since for that there are loads of formulas that I actually understand and can use). And then I go back to lon/lat using:

val lat : Double = Math.toDegrees(y/RadiusOfEarth)
val lon : Double = Math.toDegrees( x/(RadiusOfEarth*CentralLatCos) )

However, this seems to be causing quite large errors so I'm quite sure the way I project my points is flawed.

Could you help me by suggesting either and easy algorithm to project back/forth from sphere to cartesian (that keep heading&distance) or a way to interpolate a point between two other lon/lat points at a given distance ?

Also, please do provide some source material that "proves" you solution so that I can better understand it, I asked this question somewhere else and some formulas with no context where thrown at me (which doesn't really help, I can find loads of those using the internet).

Do note, I'm trying to do this for small distance, the largest of which might be something around the lines of 2-3k meters (distance between the tow original points I know). So I really have no need for a formula that's accurate on a global scale, if it has a small error but its easy to implement I'm more than fine with it.

(Note: I'm working on a map between lat 45 - 55 and lon 7-13)


texts: Table structure for a reverse geocoding database

texts:
I should start by saying that I do not have much experience with GIS so this question may be foolish. What I do have is a database for performing a lookup from a set of coordinates to find the bounding box or LSD in which they reside. At the moment the database contains several thousand tables and we run some client side code to decide on which table to query, load that entire table back to the client and filter it client side. That seems dramatically inefficient to me as we are passing back quite a lot of data. Each table defines a large area at which to look and each record defines a smaller area. The areas are defined by 4 corner points. I suppose the original design was based on a quadtree but instead of there being four children and many levels it was decided to have only two levels and many children at each level. I am looking to perform some optimization of this structure.

The database is SQL server 2005 so there is no help available from native geospatial datatypes. What I would like to know is

Is this sort of thousand+ table database common for GIS work?
Would it not be more efficient to keep the data in a single table and take advantage of some index magic? (Is it possible to index efficiently on the decimals which define each point of the bounding box such that you can find intermediate values?)
How could this structure be improved to allow for rapid queries?


texts: ST_AsGeoJSON providing unexpected values

texts:
GIS newbie here! I have been playing around with the import I did from an OSM region (that I got from here) to postgres using the osm2pgsql conversion tool. The command I used to import was this:

osm2pgsql -c -d gis --slim --style default.style -C 10240 --flat-nodes ./bayern-latest-pgsql bayern-latest.osm.pbf

Now, everything seems to have worked and I got all the data I want, except for the geometry. When I run a query like this:

> select osm_id, ST_AsGeoJSON(way), name, from planet_osm_point;

I get back values like {"type":"Point","coordinates":[1197606.06200162,6136808.52532997]} in the st_asgeojson column. These coordinates are not what I expected (lat, lon). What am I doing wrong here?


texts:
I'm using R and am trying to create a choropleth map using GISTools. My code is like:

choropleth(us_state, us_state$smallbiz1998, shades)
plot(us_interstate, add = TRUE)

us_state is a map of US states and us_interstate is a map of the interstate highways.

I want the boundaries of US states to be in gray. How can I do this?

choropleth documentation indicates that you can specify Additional parameters to be passed on to the plot method for sp.. However, when I try choropleth(us_state, us_state$smallbiz1998, shades, col = "gray") it throws an error.

(If it were changing the interstate highways to be gray, it's easy, just change the second line to plot(us_interstate, add = TRUE, col = "gray")).


texts: GeoServer export DXF with colors and labels

texts:
I have GeoServer 2.15.1 with the DXF export plugin installed and I would like to export my layers with the colors and labels shown on the GeoServer.
When I export all layers of both line and polygon types are drawn as a line and in white
In a blog post geosolutions mentions the possibility of an advanced configuration where it is possible to choose colors and labels.
https://www.geosolutionsgroup.com/blog/developers-corner-dxf-output-format-promoted-to-official-extension-for-geoserver/
Where are the color and label options for DXF file export configured?


texts: GeoServer complex WTMS/TMS/GWC trouble: No such operation gwc null service

texts:
I trying to solve a issue from replaced data dir on GeoServer. When I choose a WMTS / TMS services in GeoServer's 2.16 GUI I see this:

and 400:null too for TMS: <ows:ExceptionText>No such operation gwc null service</ows:ExceptionText>

In the log an error begins from this:

ERROR [geowebcache.GeoWebCacheDispatcher] - Request failed
java.lang.NullPointerException
    at org.geoserver.gwc.layer.StyleParameterFilter.setLayer(StyleParameterFilter.java:188)
    at org.geoserver.gwc.layer.TileLayerInfoUtil.checkAutomaticStyles(TileLayerInfoUtil.java:162)
    at org.geoserver.gwc.layer.TileLayerInfoUtil.checkAutomaticStyles(TileLayerInfoUtil.java:140)
    at org.geoserver.gwc.layer.GeoServerTileLayer.getPublishedInfo(GeoServerTileLayer.java:358)
    at org.geoserver.gwc.layer.GeoServerTileLayer.isEnabled(GeoServerTileLayer.java:280)
...

I tried to catch broken layers or styles after Tomcat's restart, but all looks ok, except this rows:

2019-10-01 12:27:16,107 INFO [storage.DefaultStorageFinder] - *** Reverting to java.io.tmpdir null for storage. Please set GEOWEBCACHE_CACHE_DIR. ***
2019-10-01 12:27:15,963 WARN [config.XMLConfiguration] - *** GWC configuration validation error: cvc-elt.1: Cannot find the declaration of element 'gwcConfiguration'.
2019-10-01 12:27:15,963 WARN [config.XMLConfiguration] - *** Will try to use configuration anyway. Please check the order of declared elements against the schema.

How can I fix it all and make wmts/tms gwc work? I have to do something with unobtained style/layer or something more using 2 last WARNs ?


texts: Minimising number of dynamic pages to map scattered points using ArcGIS Desktop?

texts:
From time to time I have to produce mapbook to show points of interest. First step to create pages, using regular mesh:

I don't like the solution because a) there are some pages with single points (e.g. page 25) sitting on the edge and b) too many pages.

First issue is easy to fix using code, - move rectangle of page extent to center of relevant points extent:

I still don't like it, it looks very crowded because number of pages remains the same.
Remember, they all end up being actual A3 paper pages in multiple copies of report !.

So I've cooked a code that reduce number of pages. In this example from 45 to 34.

I am not sure if this the best result that can be achieved,

What is the best strategy (pseudo code, publication, Python library) to shuffle through points in order to minimise number of given size rectangles to capture all of the points?
Surely, someone discovered it in game theory, military art or fishing industry

This is update to original question:

This shows real extent and page size required:

Closer zoom showing 10 out of 164 pages:

Sample Point Feature Class

Rectangle size can change as soon as it stays within the limits, i.e. smaller is fine.


texts: Extracting only number from address in QGIS Attribute Table

texts:
I want to print only the number part of the address field. I tried to use the "field calculator" with "regexp_substr" to extract just the number (only string before the first space) but couldn't figure out the right expression. I am also thinking that it could be done with labels, by showing only the first word of the address.

So "1250 Main Street" should only show "1250".


texts: batch saving 10m DTM layers using python console

texts:
I am attempting to batch save 300 10m DTM tiles all in the same Raster type 'Float 32-bit' for using to create a VRT. At present I cannot create a VRT using the original Raster layers as some of the Tiles cannot be converted due to various reasons Such as the postive NS resolution, some saved as VRT's already and others 'Float 64-bit' meaning it will not create a full VRT of these type combinations for whatever reason.

I have been able to create and save the full 300 tiles using the following python code in the console which I have gathered from this site in another question

myDir = "C:\Users\MY_DIR\"
layers = iface.legendInterface().layers()
pipe = QgsRasterPipe()
for layer in layers:
   extent = layer.extent()
   width, height = layer.width(), layer.GetRasterBand(1)
   renderer = layer.renderer()
   provider=layer.dataProvider()
   crs = layer.crs().toWkt()
   pipe.set(provider.clone())
   pipe.set(renderer.clone())
   file_writer = QgsRasterFileWriter(myDir + layer.name() + ".tif")
   file_writer.writeRaster(pipe,
                           width,
                           height,
                           extent,
                           layer.crs())

The problem is that once created in my directory folder, the .tif's created have got 4 raster bands and the heights of my Terrain are all lower than what the original raster is. Is there something I could add to this code that would allow me to keep the original Band one values of the raster corresponding to the true height of the terrain?


texts:
I am trying to use ArcGIS Pro to create a 3D scene using an elevation surface (.tif) generated from drone imagery. The tif is projected and contains the elevation information, but I cannot get ArcGIS pro to accept the layer in 3D - only in 2D.
I can use the raster as an elevation surface in the scene:
Right-click Scene in Contents Pane > Properties > Elevation Surface > ground > navigate to raster layer
But, the layer doesn't remain 3D when I publish the Web Scene to ArcGIS Online.
How do I use the raster as the actual 3D layer?


texts: pgr_pointsaspolygon returns a self intersecting polygon which does not cover underlying points

texts:
I am using pgrouting 2.0 library to create an alpha shape around my set of points.
It works nearly perfect but I have a case where pgr_pointsaspolygon returns a self intersecting polygon which does not really cover the underlying points in the table:

This is how i call the pgr_pointsaspolygon function :

SELECT 1 as a,* FROM  pgr_pointsaspolygon('SELECT id, ST_X(the_geom) AS x, ST_Y(the_geom) AS y FROM test')

Is there anyway i can correct the result?

A bug was report here: https://github.com/pgRouting/pgrouting/issues/235


texts: Adding feature class and applying symbology in ArcGIS Pro using ArcPy?

texts:
Using an ArcPy script run from a toolbox, I want to import a feature class from a geodatabase into an ArcGIS Pro 2.1.2 map and then symbolise that layer based on a lyrx file stored on disk.
Here is a simplified version:
import arcpy, os
aprx = arcpy.mp.ArcGISProject('CURRENT')
m = aprx.listMaps('*')[0]
outputWorkspace = r"C:\Path\To\MyProject.gdb"
arcpy.env.workspace = outputWorkspace
out_fc = "fc_data"
out_fc_p = os.path.join(outputWorkspace, out_fc)
ref_lyrx = r"C:\Path\To\symbol_ref.lyrx"
m.addDataFromPath(out_fc_p)
out_fc_lyr = m.listLayers()[0]
arcpy.ApplySymbologyFromLayer_management(out_fc_lyr, ref_lyrx)

The script completes without error but no visible change is made to the symbology.
Is this a bug?
I've reviewed many similar questions but have been unable to resolve this.
I tried also first converting the fc to a layer using arcpy.MakeFeatureLayer_management().getOutput() before adding to the map with m.addLayer() and then using arcpy.ApplySymbologyFromLayer_management() but get exactly the same outcome.