texts: Are Xbase-dbf-files still common standard for shapefiles?

texts:
For a Java application and for my personal understanding, I am exploring how to read and write shapefiles (at least shx, shp, dbf). The main purpose is to keep the resulting Shapefile-class and it's complimentary classes relatively lightweighted but with a backdoor to be able to react on changes of the different file formats as well.

The spatial information is well documented and clearly structured. This goes for the data table (dbf) to some extend as well. However, I am struggling with the different versions of dbf-files at the moment. I know that my version of QGIS Lisboa saves dbf-files as Xbase-files in a fashion which is described here. At the current stage, my Java classes are able to handle those. On the other hand I have read that different dBase-formats might be in use too. The respective file-specification is documented on e.g. this website.

Is the Xbase-format still dominant in comparison to the more complex but also more flexible dBase-formats?


texts:
I don't find a way to see the 11 points of my set of data that are just meters from each other as separate points in GeoServer and I don't know why. I would have assumed that GeoServer finds the correct resolution through calculating the boundary box but that doesn't seem to help. Maybe you can reproduce it given the following information:
I have the following table in PostgreSQL 13.2 with PostGIS 3.1 (currently on windows 10 but setting up a linux server to move to that)
                          Tabelle ╗public.samplepoints½
   Spalte    |         Typ          | Sortierfolge | NULL erlaubt? | Vorgabewert
-------------+----------------------+--------------+---------------+-------------
 sampleid    | integer              |              | not null      |
 sampletag   | smallint             |              |               |
 sampletype  | smallint             |              |               |
 samplegeom  | geometry(Point,3035) |              |               |
Indexe:
    "treesamples_pkey" PRIMARY KEY, btree (sampleid)

with the following data:
 sampleid | sampletag | sampletype |                     samplegeom
----------+-----------+------------+----------------------------------------------------
        1 |     20523 |        949 | 0101000020DB0B0000EF8DAD60850123409ECB585D1C7A4940
        8 |     20530 |        949 | 0101000020DB0B00008DC19DA4AD0123407C4B488A197A4940
        2 |     20524 |        837 | 0101000020DB0B000057BB699C8D012340FC7DBFCB1E7A4940
        3 |     20525 |        949 | 0101000020DB0B0000E86C4F7985012340F84EED881A7A4940
        4 |     20526 |        188 | 0101000020DB0B0000825218C28F012340F211BC171B7A4940
        5 |     20527 |        356 | 0101000020DB0B0000EC2ED81E99012340BF414BBC1B7A4940
        6 |     20528 |        578 | 0101000020DB0B0000F8CFA08F9C0123406BF6A09D1C7A4940
        7 |     20529 |        554 | 0101000020DB0B00005B4CFA9FA60123400D74F77D1B7A4940
        9 |     20531 |        578 | 0101000020DB0B00007574EB5CA90123401EED730D177A4940
       10 |     20532 |        239 | 0101000020DB0B000072F269319B012340B350770D167A4940
       11 |     20533 |        748 | 0101000020DB0B0000B38462578D012340EB350836187A4940

Then I create a workspace "sampleworkspace" in GeoServer and a store "samplestore" with access to PostGIS.
Finally, I create a layer from sampleworkspace:samplestore in which I publish samplepoints from PostGIS.
I keep the referencing system as EPSG 3035 and let the layer's bounding boxes be calculated from the data which results in
Min X        | Min Y        | Max X        | Max Y
9,5029707157 | 50,953797992 | 9,5032779162 | 50,954064816
Lat/Lon
Min X              | Min Y              | Max X              | Max Y
-29,08689230212812 | 12,994046968698024 | -29,08689229881456 | 12,994046972307304

So my 11 points are very close to each other.
A preview of the layer as OpenLayers only shows me one point and the scale is set to 1:2 which means I cannot zoom in any closer, only zoom out further.
When I click on that one point, it gives me all 11 points:
samplepoints
fid             | sampletag | sampletype
samplepoints.1  | 20523     | 949
samplepoints.8  | 20530     | 949
samplepoints.2  | 20524     | 837
samplepoints.3  | 20525     | 949
samplepoints.4  | 20526     | 188
samplepoints.5  | 20527     | 356
samplepoints.6  | 20528     | 578
samplepoints.7  | 20529     | 554
samplepoints.9  | 20531     | 578
samplepoints.10 | 20532     | 239
samplepoints.11 | 20533     | 748

Is there any way to see the 11 points as separated points, so to go "beyond 1:2" or to adjust the scaling in a different way?


texts: Seeking host for free (little) wms server (for fantasy project)

texts:
I am programming a website (Node server) for a friend of mine who created a fantasy world, to promote his books and games.

As a former GIS tech, I really would like to setup a map of his world (using his drawings). I don't have any worries about projection since I'll be using Mercator, and nor about rendering since I can choose between OpenLayers and Leaflet.

My problem is: I don't own a server, nor any computer that could be used as, I rely on Heroku for hosting my apps and websites freely.

Does anyone know of a solution, some kind of WMS host allowing to upload at least georeferenced rasters?

I considered installing Geoserver or Mapserver but I really don't think free plans like I have on Heroku or any other webhost would allow me to install such softwares.

I considered also using npm packages directly in my website backend, but all I can find seems to be vector and no raster

There are about two dozen standard paper (A4) drawings to upload, I am open to any format.

If I find nothing else, I will vectorise the drawings and use a node server.


texts: How to enable adding both feature classes / tables and layers from TOC to a script tool parameter input?

texts:
I have developed a script tool which will allow a user to specify input feature class or table and then will process it.

Surprised I haven't noticed this before: I cannot find which data type my input features parameter should be to let users specify both feature class / tables (by drag-n-drop from Catalog window or browsing to a place on the disk) and layers from the table of contents (TOC) in case they have already added a feature class / table to the map document.

Dataset and Feature Class type will let me to browse to tables/feature classes and feature classes only, respectively. Layer will let me specify the added feature layer from the TOC only.

Core GP tools such as Frequency and many others, in contrast, will let you drag-n-drop feature classes/tables and choose a layer/table in case you have some added to the TOC.

I've read through the list of GP data types from the Esri Help, but cannot find any that will help me. Am I hitting one of the limitations of the custom script tools? I am on 10.2.1.


texts: Old MIller projections Bipolar Oblique Conic Conformal/Oblated Stereographic definitions?

texts:
I'm attempting to georeference and warp some old toposheets that are various old projections including the following:

Miller Oblated Stereographic Projection
Bipolar Oblique Conic Conformal

What's the best approach to acquiring the definition for this projection?  Or do I just have to use ArcGIS which seems to have these defined, note I have ArcGis available, but not for a while.


texts:
I have a csv with a list of the world sea ports and their coordinates.

The coordinates are in this format:

5234N 00226W

How can i convert this to decimal Latitude / Longitude ?

Data source


texts: What determines whether you can open an attribute table inqGIS?. I'd like to edit an attribute table for a raster in qGIS, but it won't me open the table. There's plenty in the help about editing attribute tables, but it doesn't say when you can/can't open them


texts:
I created a Python plugin in QGIS & when i am running my plugin it printing the data on python console. But i want that when i run my plugin it will print the data on QGIS map canvas. How can i do this?
Here is my run method which is going to print lat/long per second on Map Canvas. you just tell me how can i get lat/long on Map Canvas instead of printing it in python console.
Remember Python version 3.6 QGIS version 3.2 & this code working fine no error

def run(self):
    print('hello')
    """Run method that performs all the real work"""

    PORT_NUMBER = 5000
    SIZE = 1024

    hostName = gethostbyname('192.168.225.140')

    mySocket = socket(AF_INET, SOCK_DGRAM)
    mySocket.bind((hostName, PORT_NUMBER))

    print ("Test server listening on port {0}\n".format(PORT_NUMBER))
    print ("Hi i am linux computer")

    while True:
        (data, addr) = mySocket.recvfrom(SIZE)
        print (data)

    self.dlg.show()
    # Run the dialog event loop
    result = self.dlg.exec_()
    # See if OK was pressed

    if result:
        # Do something useful here - delete the line containing pass and
        # substitute with your code.
        pass


texts:
I have a buffer and need to filter a multipoint layer (QGSVectorLayer) so that I would get only features intersecting the buffer. What I've found:

a) setFilterRect()

This doesn't accept polygon arguments, just rectangles. I've found no polygon alternative. I don't want to decrease accuracy by converting my buffer to its MBR.

b) setSpatialFilter()

I must be using it wrong way, or perhaps it's a part of an old API which is no longer used. In the example it's applied directly on the layer, but while trying this, I get an error:

AttributeError: 'QgsVectorLayer' object has no attribute 'SetSpatialFilter'

c) select()

I've found that such a function exists, but my google-fu is too weak to find detailed documentation on it or whether it's suitable for my needs. This link seems to say that it's deprecated/abandoned now, though I might have misunderstood it.

There already are few similar questions, even almost duplicates. One is basically "I did this in ArcPy, is it possible in PyQGIS?", and the only answer (in comments) is "it is possible" with no suggestion how to do it. The others provide some workarounds, either through looping and calling intersects() or through QgsOverlayAnalyzer. Neither of these two seem to be a clean solution, so even though I will try the intersects() way, I don't consider it a solution unless someone proves that there's no better way.


texts: Convert tabulate area's output format into a text file (.txt or .csv) in the ModelBuilder in ArcGIS?

texts:
When I use the Tabulate Area tool, the original output's format is .dbf.

What can I do to export the output as a text file (.txt or .csv) in the ArcGIS ModelBuilder environment?


texts:
I am trying to add new field to my table based on two other existing fields.
In the field Calculator it's quite simple:

When I am trying to run as a Python code I got into trouble, Here's my code:

vLayer = iface.activeLayer()
vLayer.dataProvider().addAttributes([QgsField('1m_Cases', QVariant.Int)])
vLayer.updateFields()
for f in vLayer.getFeatures():
        f.setAttribute(f.fieldNameIndex('1m_Cases'), f(['J_TotalCases'] *1000000)/ f['Pop'])
vLayer.commitChanges()


texts: Turn line labels in QGIS to a new point shape layer

texts:
I have a set of contour lines in shape format that I have labeled in QGIS. Now I'd like to export these labels as point features with attribute data containing the text of the label, preferably in shape format. Is this possible?

Also, if there is another open source tool to do the labelling which makes it easier to create a point layer from I could use that one.

(The points will end up in MAPublisher / Illustrator)


texts: Export GPS file in Lambert 93 for Illustrator manipulation

texts:
In Illustrator, I use some French IGN maps (1/50.000 scale) Lambert 93.
I want to put on some GPX files to the right position for sure.
I use QGIS to import the GPX file, I modify the reference to Lambert 93 (EPSG:2154). I export it in DXF file. When I open it in Illustrator, it's impossible to have the GPX file in the right place on my map. Any idea ?


texts: Creating ArcGIS toolbar with new tools using Python?

texts:
I have written a number of fairly-simple python functions using arcpy which help a user significantly when performing a certain workflow using various tools (both ArcGIS toolbox tools and external software).

I'd really like to bind some of these functions to buttons on a new ArcGIS toolbar, as they are not really geoprocessing tools, but more like 'Workflow Management Tools', as they do things like show and hide various layers, change symbologies etc, and generally make it easier for the user to perform specific tasks.

Is there a way to do this? I have looked in the documentation and have found various things about how to create toolbar buttons using ArcObjects, and I have found this question about accessing ArcObjects from Python, but I'm not sure how easily these can be combined.

Does anyone know if this can be done, and if so, how?


texts: Merging two spatial objects in R: spCbind gives error "row names not identical"

texts:
I have two SpatialPolygonsDataFrames. One is the LSOA shapefile of Great Britain (i.e., Great Britain divided into ~41,000 regions with approximately equal population densities), the other is a bunch of regions in Great Britain (coal fields). I want to add the coal field info to the LSOA object, in order to have a variable added to the LSOA object that tells me whether an LSOA region falls into a coal field or not.

I got the shapefiles from here:

LSOA regions:
https://www.statistics.digitalresources.jisc.ac.uk/dataset/2011-census-geography-boundaries-lower-layer-super-output-areas-and-data-zones/resource-1
Coal mines: https://data.gov.uk/dataset/coal-mining-reporting-area

This is how I attempted to do the merge:

# read in LSOA regions
LSOA <- readOGR(dsn = "Infuse", layer = "infuse_lsoa_lyr_2011")
# transform map to OSGB projection
bng = '+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +datum=OSGB36 +units=m +no_defs'
LSOA.bng <- spTransform(LSOA,bng)

# read in coalmining regions:
coalmine <- readOGR(dsn = "coalmining_shapefiles", layer = "Coalfield_polygon")
# transform map to OSGB projection
coalmine.bng <- spTransform(coalmine,bng)

# merge data
coalmines.LSOA.over <- over(LSOA.bng, coalmine.bng)
LSOA.coalmines <- spCbind(LSOA.bng,coalmines.LSOA.over)

This is the error message I get after the spCbind command:

Error in spCbind(LSOA.bng, coalmines.LSOA.over) : row names not identical

What did I do wrong?


texts: GDAL Coordinate Transformation: NoneType Error

texts:
I am trying to transform the coordinates of the shapefile and have been getting mixed results from this script when used for different shapefiles.
from osgeo import ogr, osr
import os

driver = ogr.GetDriverByName('ESRI Shapefile')

# input SpatialReference
inSpatialRef = osr.SpatialReference()
inSpatialRef.ImportFromEPSG(28356)

# output SpatialReference
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromEPSG(4326)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# get the input layer
inDataSet = driver.Open(r'Rail_network.shp')
inLayer = inDataSet.GetLayer()

# create the output layer
outputShapefile = r'Rail_network_WGS84.shp'
if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)
outDataSet = driver.CreateDataSource(outputShapefile)
#for point replace Polygon with Point or MultiPoint, for line replace with LineString or MultiLineString
outLayer = outDataSet.CreateLayer("basemap_4326", geom_type=ogr.wkbMultiLineString)

# add fields
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
inFeature = inLayer.GetNextFeature()
while inFeature:
    # get the input geometry
    geom = inFeature.GetGeometryRef()
    # reproject the geometry
    geom.Transform(coordTrans)
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    outLayer.CreateFeature(outFeature)
    # dereference the features and get the next input feature
    outFeature = None
    inFeature = inLayer.GetNextFeature()

# Save and close the shapefiles
inDataSet = None
outDataSet = None

#create prj file
spatialRef = osr.SpatialReference()
spatialRef.ImportFromEPSG(4326)

spatialRef.MorphToESRI()
file = open('Rail_network_WGS84.prj', 'w')
file.write(spatialRef.ExportToWkt())
file.close()

The error was:
line 44, in geom.Transform(coordTrans) AttributeError: 'NoneType' object has no attribute 'Transform'`

I tested the shapefiles on QGIS and they were showing up fine there. Can anyone identify what has gone wrong?


texts: How to create contour lines from TIFF (raster with more than 100 bands)

texts: I need to do batch processing from a raster with 120 bands (GeoTIFF). I need 120 contour line shapefiles, but I get only one. I had no problem when the result is another raster file with all the bands, but the contour line is a shapefile, I don´t know if this is the problem...


texts: Setting color based in quotient and remainder of data with Geoserver WMS SLD?

texts:
I need to style a raster WMS layer based on a math operation over the raster data. In this case I need to set for example one output channel as a linear ramp from the quotient of the raster data by a constant, and another channel from its remainder.

Is this possible using the SLD file?


texts: How to automate an insert in a PostgreSQL database using Python

texts:
I'm developing a plugin in qgis that makes a selection, identifies the IDs of the selected elements, and wanted it to insert the selected IDs into the PostgreSQL DB.

I have this code to do an insert into the PostgreSQL database from python:

import psycopg2
conn = psycopg2.connect("\
        dbname='postgres'\
        user='postgres'\
        host=''\
        port=''\
        password=''\
");
c = conn.cursor()
c.execute("INSERT INTO (ID) VALUES ('')")
c.execute("SELECT * FROM ")
conn.commit

records2 = c.fetchall()
import pprint
pprint.pprint(records2)

But i wanted the insert to be automated, for exemple, it is not necessary for me to write which IDs are selected,
but be the plugin to identify and insert the IDs based on the selection i made.


texts: Labeling in TileMill with CartoCSS - Add line break between variables

texts:
I am labeling some points in TileMill using CartoCSS.

The labels are a composite of two variables

text-name: [NAME]+' '+[TYPE];

The rule above inserts a space between the two variables. How can I add a line break?

I have tried

text-name: [NAME]+
+[TYPE];

.

text-name: [NAME]+'
'+[TYPE];

.

text-name: [NAME]+'_'+[TYPE];
text-wrap-character: '_';

All to no avail. Any suggestions?


texts: Automatically delete services after a set period of time

texts:
Is it possible to automatically delete services after a set period of time? Is the username of the publishing user and the date the service was published saved anywhere? We are trying to implement a system to check for services that have been on the server for three months, identify the user that published them, and automatically email the user to find out if the service should be saved any longer. If yes, remind the user in another month. If no, delete the service. If no response, delete the service.

Any ideas for how to approach implementing this? Just looking for overall strategies here. Thanks


texts: Geopdf map not displaying coordinates in Adobe Acrobat Geospatial Location Tool.

texts: I georeferenced a map using a data layer from ArcGis online as a base map. I transformed the now georeferenced map to the WGS 1984 World Mercator projected coordinate system. Lastly, I exported the map from ArcGIS as a PDF. I exported all data layers. However when I open in Adobe and use the Geospatial Location tool the lat and long bar remains blank. What am I doing wrong?


texts: Calling a QgsProcessingAlgorithm from a Processing python script in QGIS 3.0?

texts:
I am writing a Processing python script based on the template provided here.

I am creating lines (destined to the 'OUTPUT') and I would like to clip them according to a polygon entered as the 'INPUT' parameter).
I can't figure out how to call the Algorithm. For what it's worth, this doesn't work:
QgsProcessingAlgorithm('qgis:clip').run({'INPUT': self.OUTPUT, 'OVERLAY': self.INPUT})

I have tried several ways but couldn't get far as I get the error:
qgis._core.QgsProcessingAlgorithm represents a C++ abstract class and cannot be instantiated


texts:
With QGIS 2.8 I need to create some composers based on a template, and save them using Python.

I've found how to load a template in a composition, but I can't find how to save this composition in the list of composers, nor how to create a QgsComposerView.

Do you know how to do so, or tell me where to look?


texts: Using WofE (weight of evidence) in ArcMap

texts:
I need to use "weight of evidence" in a project about floodrisk
I heard that this model can be implemented in ArcMap but I did not find any videos in YouTube.
Can anyone help me?


texts: Is there a way to create multiple points in one layer of a Cartodb map and edit them each separately?

texts: I am making a Cartodb map that is supposed to show different pictures of a park in popups at points along a park path.  I started by making a different layer for each point but by the time I created 4 layers it looked like it was going to stop allowing me to create any more.  I have around 30 unique popups that I need to make.  Is there a way to create multiple points in one layer and edit them each separately?


texts: How to call processing toolbox in Qgis' python console

texts:
I have been trying to use processing toolbox via Python console in QGIS 2.0.1. I am able to write the following commands:

I am particularly interested to run 'r.statistics' from GRASS module.

My question is how to arrange the input parameters sequentially and what does each parameter mean? For example, in the above screenshot, when I type the help, I can see the parameter inputs,
processing.runalg('grass:r.statistics', base, cover, method, -c, GRASS region parameter, Grass region cell size, output)

I can understand that 'base = input raster', 'method = avg, mean, std' and 'output=outputfile to eb declared'.
But I am not sure what are 'cover=', '-c','GRASS region means'? How can I get further help on that? If I were to run the function 'r.statistics' via processing toolbox, how should the exact/correct command look like?
processing.runalg('grass:r.statistics', 'C:\mydata\dem.tif', cover, '1', -c, GRASS region parameter, Grass region cell size, 'C:\mydata\outputdem.tif')

Can anyone help me what to put in place of rest of the parameters?
Thanks in advance.


texts: Running ArcGIS Desktop in virtual machine?

texts:
I'm trying to run ArcGIS Desktop 10.0 in a Windows 7 (32bit) virtual machine. I don't have any native Microsoft Windows available. I'm using VirtualBox from Oracle for virtualization in a native linux environment (Sabayon).

I'm having exactly this issue: ArcGIS Explorer Map Display Issue

The proposed solution is to update video card driver. But in VirtualBox there are driver called "Guest additions" and they are already up-to-date.

I know some of my fellow students are working with ArcGIS in VMs and they seem not to have any kind of issues. What could it be?

3D- and 2D-Video-acceleration are enabled.

It's a native Linux machine, running a virtual Windows 7 machine - 32bit.
I'm aware of ESRI's system requirements. I allocated 256MB video RAM and enabled 2D and 3D hardware acceleration - as seen in the screenshot. Anyways, I've been running this analyzation of my system requirements from inside the machine, this is the result:

Now I got to figure out how to fix it. Video memory is okay, but the rest I've never heard of...

I installed the experiemental WDDM drivers as @scw suggested, now it looks like that:

I have all recommended attributes now, but it still says:

  Unfortunately, your Video Card does not meet this requirement.

In addition, ArcGIS crashs now after adding any layers.


texts: Cartodb TypeError: layer.getSubLayer is not a function

texts:
Can anyone tell me why I am getting this error? It happens when I try and filter by a layer action. I also tried to do the cartodb tutorial on "Toggle the Map View" following it line for line and am getting the same error.

var map;
function init(){
  // initiate leaflet map
  map = new L.Map('map', {
    center: [39,-95],
    zoom: 4,
    minZoom: 3
  })

  L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', {
      attribution: '&copy; OpenStreetMap contributors'
  }).addTo(map);

  var layerUrl = 'https://myusername.cartodb.com/api/v2/viz/88afc111-16cf-11e5-a0d9-0e9d021ea90d/viz.json';

  var layerOptions = {
            query: "SELECT * FROM mytable_1",
            legends:false,
            cartodb_logo: false
  }

  cartodb.createLayer(map, layerUrl, layerOptions)
   .on('done', function(layer) {

    map.addLayer(layer);
  }).on('error', function() {
    //log the error
  });

  var LayerActions = {
    all: function(){
        layer.getSubLayer(1).setQuery("SELECT * FROM mytable_1");
        return true;
      },
    one: function(){
        layer.getSubLayer(1).setQuery("SELECT * FROM mytable_1 WHERE overall_rating ILIKE '%1%'");
        return true;
      },
    two: function(){
        layer.getSubLayer(0).setQuery("SELECT * FROM mytable_1 WHERE overall_rating ILIKE '%2%'");
        return true;
      },

    three: function(){
        layer.getSubLayer(0).setQuery("SELECT * FROM mytable_1 WHERE overall_rating ILIKE '%3%'");
        return true;
      },
    four: function(){
        layer.getSubLayer(0).setQuery("SELECT * FROM mytable_1 WHERE overall_rating ILIKE '%4%'");
        return true;
      },
    five: function(){
        layer.getSubLayer(0).setQuery("SELECT * FROM mytable_1 WHERE overall_rating ILIKE '%5%'");
        return true;
      },
    na: function(){
        layer.getSubLayer(0).setQuery("SELECT * FROM mytable_1 WHERE overall_rating ILIKE '%NA%'");
        return true;
      }
  }
  $('li').click(function(){
    $('li').removeClass('selected'); $(this).addClass('selected');
    LayerActions[$(this).attr('id')]();

  })
}


texts:
I'm trying to use some of the GRASS tools from the processing toolbox in QGIS. I get the following Python error:

An error has occured while executing Python code:

Traceback (most recent call last):
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\gui\ProcessingToolbox.py", line 181, in executeAlgorithm
    message = alg.checkBeforeOpeningParametersDialog()
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Algorithm.py", line 511, in checkBeforeOpeningParametersDialog
    msg = Grass7Utils.checkGrass7IsInstalled()
  File "C:/PROGRA~1/QGISCH~1/apps/qgis/./python/plugins\processing\algs\grass7\Grass7Utils.py", line 361, in checkGrass7IsInstalled
    cmdpath = os.path.join(path, 'bin', 'r.out.gdal.exe')
  File "C:\PROGRA~1\QGISCH~1\apps\Python27\lib\ntpath.py", line 96, in join
    assert len(path) > 0
TypeError: object of type 'NoneType' has no len()

Python version:
2.7.4

QGIS version:
2.4.0-Chugiak

Any ideas on how I can fix this? I've tried a re-installing but haven't had any luck - same error keeps on reappearing.


texts: Using OGC Filter with GeoNetwork to make CSW request?

texts:
I have a GeoNetwork catalog with a lot of ISO 19115 metadata. I want to use ogc filters to search into my catalog through CSW services (GetRecords).

So, I use http POST requests. The post data respect GeoNetwork XML schema.

URL :

http://localhost:8080/geonetwork/srv/eng/csw

PostData:

<?xml version="1.0" encoding="UTF-8"?>
<csw:GetRecords service="CSW" version="2.0.2" xmlns:csw="http://www/opengis.net/ogc"
                                              xmlns:ows="http://www/opengis.net/ows"
                                              xmlns:gml="http://www/opengis.net/gml">
  <csw:Query typeNames="csw:Record>
    <ElementSetName typeNames="csw:Record>full</ElementSetName>
    <csw:Constaint version="1.1.0">
      <ogc:Filter>
        <ogc:Intersects>
          <ogc:PropertyName>ows:BoundingBox</ogc:PropertyName>
          <gml:Envelope>
            <gml:lowerCorner>5 42</gml:lowerCorner>
            <gml:upperCorner>8 50</gml:upperCorner>
          </gml:Envelope>
        </ogc:Intersects>
      </ogc:Filter>
    </csw:Constaint>
  </csw:Query>
</csw:GetRecords>

It works for logical requests (title is ...) but I do not have any result for spatial requests. The response status code is 200 but I have 0 record returned even if I have metadata in the bounding box (Envelope).

Is my post data correct?

Do the lower corner and upper corner correspond to "long lat"?

Do you know good documentation about ogc filter?

I also use ogc filter to update metadata in my GeoNetwork catalog (transaction CSW services) but I can filter only by some tags.

Do you know how to add filters to GeoNetwork? Have I to modify the config-cw.xml GeoNetwork file?


texts: For the purposes of spatial interpolation, I have created a grid using "Create Fishnet" in ArcMap10.2 (I use the "label points", which have a spacing of 0,00833333 degrees). For every point (label point) I want to have the mean elevation (aspect,curvature) in 800m. around the point. Any ideas how this could be performed?


texts: Taking attribute from multiple points within fixed distance

texts:
I have a list of point features each with a unique "ID". I need to run a fixed distance buffer on all 70,000 points to establish which points are within a fixed distance of each other than return the unique "ID"s for each point. These need to be within columns in the attributes table.

So point 1 (id-A) might have point 9 (id-x) and Point 87 (id-c) within 1 km. I need the attributes table of point 1 to list id-A, then id-x then id c, all in columns for export to Excel where I can concentrate them with comers in between.


texts: Preparing raster Image for correlation in ENVI 5.5

texts: I am working with two raster image of the same location with the same extent, both projected to the same coordinate system and having the same cell size. I want to correlate the raster value on the images using ENVI 5.5. I have classified the raster value of these images into two classes but I am interested in the correlation of only one class on both images. How do I extract the only class of interest on both images, thereby leaving other class as area of no data, before running the correlation?


texts: QGIS 2.18.9: buffers too large in WGS84/UTM zone 34S

texts:
my original vector layer was saved in WGS84. In order to be able to calculate distances in meters I have saved the layer in WGS84/34S. I thought that the unit of measurement would be meters now. But when I create 50 meter buffers around buildings they come out way too large. When measuring distances they are given in milimeters instead of meters.

Which steps should be taken in order to get to a unit of measurement of meters?


texts: Problem with intersection between two Geopandas Dataframes

texts:
I have 2 shape files one "country.shp" (containing the tiles geometry covering the whole country) & other "district.shp"(containing one particular district). Now I want to get the tiles information getting intersected with the district polygon (similar to search by location in ArcGIS) and the final intersected geodataframe should contain the attributes of "country.shp" (containing only intersected rows from "country.shp"). Up till now I have used the following code but it is working fine only for some district shape files (projections are same for both the shapefiles). Can someone help me what is wrong with the following code or suggest some other alternatives?

import geopandas as gpd

f1 = "/home/geo/country.shp"
f2 = "/home/geo/district.shp"
data1 = gpd.read_file(f1)
data2 = gpd.read_file(f2)

import shapely.speedups
shapely.speedups.enable()

pip_mask = data1.within(data2.loc[0, 'geometry'])

pip_data= data1.loc[pip_mask]


texts:
I'm facing problems exporting a pointcloud in geo-coordinates in my software. I'm using liblas to export the point data.

I have GeoData to display but upon loading it I subtract an offset to keepo the values low.
Initially when exporting I set the scale and offset in the liblas::Header to 1e-4 scale and 0 offset (for testing only). This gives me a decent, quite dense point cloud (see first screenshot below).
Now I'm experimenting with exporting to the global coordinate system, so I set the scale to 1 (my coordinates are in meters) and the offset to the offset with respect to the UTM coordinate system, so some value like (354000,5403000,400).

In CloudCompare the right sale/offset is shown and the point cloud looks reasonable. However, It's displayed really sparse (second screenshot). Any ideas where this might come from? Could it be an issue when exporting or a loss of exactness in CloudCompare? The number of points / point cloud size in MB is exactly the same.


texts:
I am looking at collecting spatial information from a range of stakeholders to identify opportunities within a landscape.  However, the data that we are collecting are not presently mapped and so we wish to use an interactive method to get this information onto a map.

Which online websites/methods would people recommend for circulating to a non-technical audience to try and get them to map what they know?

Where the emphasis is on ease of use by the stakeholder.


texts: Why isn't StackBuilder installation of PostGIS working?

texts:
I've installed PostGIS (ver. 2.0.1) on my Mac using EnterpriseDB's StackBuilder, which I also used for my PostgreSQL (9.2.2) installation. Although StackBuilder says that PostGIS is now installed, and although I have a PostGIS directory in /usr/local/PostgreSQL/9.2/, it doesn't appear to be actually installed (see queries below).

gregory=# select postgis_full_version();
ERROR:  function postgis_full_version() does not exist
LINE 1: select postgis_full_version();
               ^
HINT:  No function matches the given name and argument types. You might need to add explicit type casts.

gregory=# create extension postgis;
ERROR:  could not open extension control file "/usr/local/PostgreSQL/9.2/share/postgresql/extension/postgis.control": No such file or directory

and select * from pg_available_extensions; does not return a record for PostGIS.

Does anyone know what's up with this and how to get it to work?


texts: Crawling directory and listing all datasources using ArcPy?

texts:
I'm editing a script to try and create a list of all datasources within MXD's. I've got the script to crawl the file tree and find all the info, but I would like it to output the data for each layer/datasource on a single line, so it is easy to sort in excel.

Can someone take a look at the code I have so far and let me know how to accomplish that?

import arcpy, datetime, os

try:

arcpy.gp.overwriteOutput = True

#Read input parameters from GP dialog
folderPath = arcpy.GetParameterAsText(0)
output = arcpy.GetParameterAsText(1)

#Create an output file
outFile = open(output, "w")

#Loop through each MXD file
count = 0
for root, dirs, files in os.walk(folderPath):
    for filename in os.listdir(folderPath):
      fullpath = os.path.join(folderPath, filename)
      if os.path.isfile(fullpath):
        if filename.lower().endswith(".mxd"):

            #Reference MXD
            mxd = arcpy.mapping.MapDocument(fullpath)
            count = 1

            #Write MXD data to file
            MapDoc = os.path.basename(mxd.filePath)
            MapDocPath = mxd.filePath

            #Reference each data frame and report data
            DFList = arcpy.mapping.ListDataFrames(mxd)
            for df in DFList:
                #Format output values
                if df.description == "": descValue = "None"
                else: descValue = df.description

                #Reference each layer in a data frame
                lyrList = arcpy.mapping.ListLayers(mxd, "", df)
                for lyr in lyrList:
                    lyrName = lyr.name
                    if lyr.supports("dataSource"):
                        lyrDatasource = lyr.dataSource
                    else: lyrDatasource = "N/A"

            outFile.write("Map Document, MXD Path, DataFrame List, Description, Layer name, Layer Datasource")
            str = ",";
            seq = (MapDoc, MapDocPath, DFList, descValue, lyrName, lyrDatasource);
            outFile.write(str.join(seq))

            del mxd

outFile.close()

#Open resulting text file
os.startfile(output)

#Delete variables that reference data on disk
del folderPath, output, outFile, fullpath

except Exception, e:
import traceback
map(arcpy.AddError, traceback.format_exc().split("\n"))
arcpy.AddError(str(e))

The error I recieve when I try to run this is as follows:

\MultiMXDReport_crawl.py", line 60, in <module>
outFile.write(str.join(seq))
TypeError: sequence item 2: expected string or Unicode, list found
<type 'exceptions.TypeError'>: 'str' object is not callable
Failed to execute (MultiMXDCrawler).


texts: How can I finish my incomplete GPX tracks?

texts:
I have some incomplete GPS recordings. For whatever reason, the tracks weren't recorded , so there are some gaps in it. I could easily use the satellite images using QGIS software.

So I would like know how can I get them in order to alter an recorded GPX tracks without losing the timestamps.
I want to:

add tracks (track points)
split tracks
merge tracks
add timestamps to manual added tracks/ track points based on the original ones
snaps tracks / track points to features (mainly roads) in osm or google maps.
export the altered records so that geosetter can read them flawless.


texts: How to set band names of multi-band GeoTIFF asset upon upload to Google Earth Engine

texts:
I have created a multi-bands GeoTIFF using xarray in Python following this answer. When I display this raster or open it in QGIS, all the bands are named following the variable they represent:

However, when I upload the file as a Google Earth Engine raster asset, all band names are automatically renamed to b[1-9]:

Is there a way to preserve the band names when importing a GeoTIFF? Or am I missing something in the way multi-bands GeoTIFFs are defined?
Here is the code I use. I first extract data from a NetCDF file and then write a xarray dataset:
import xarray as xr
import rioxarray as rio
import pandas as pd

#%%
rootDir = '/Users/seb/Documents/WORK/Projects/Vegetation/Fall3D/CC2011/'
dataVar = 'max'

ncfile = xr.open_dataset('{}CC2011_{}.res.nc'.format(rootDir, dataVar))

className = ['lapilli-01', 'lapilli-02', 'lapilli-03',  'coarse_ash-01',  'coarse_ash-02', 'fine_ash-01','fine_ash-02','fine_ash-03',    'fine_ash-04',    'fine_ash-05',   'aggregate-01']

data = ncfile['tephra_grn_load']
# Get last date
data = data[-1,:,:]
#Define lat/long and remove the time dimension
data = data.rio.set_spatial_dims('lon', 'lat').drop('time')
#Check for the CRS
data.rio.set_crs("epsg:4326")

# Define dataset
xr_dataset = xr.Dataset()
xr_dataset['total'] = data

for bin in range(0, len(className)):
    dataBin = ncfile['tephra_grn_load_bin']
    # Get last date
    dataBin = dataBin[-1,bin,:,:]
    #Define lat/long
    dataBin = dataBin.rio.set_spatial_dims('lon', 'lat').drop('time').drop('tephra_bin')
    #Check for the CRS
    dataBin.rio.set_crs("epsg:4326")
    # Set name
    dataBin.name = className[bin]
    # Add to dataset
    xr_dataset[className[bin]] = dataBin

ncfile.close()

# Save
xr_dataset.rio.set_crs("epsg:4326")
xr_dataset.rio.to_raster('{}CC2011_{}_load.tif'.format(rootDir, dataVar))


texts: Distinction between Difference and Clip tools in QGIS?

texts: What is the primary distinction between the Difference and Clip geoprocessing tools in QGIS?


texts: Java GIS framework with custom shaped markers

texts: I am completely new to GIS theme and I found this board so I hope I can get good advices from experienced users. I need to find an open-source Java framework/library/api that would generate an map image with multiple coordinates markers. In particular I need to customize the shape of the markers (i.e. making the marker circle size bigger or smaller depending on external additional data). Moreover, it would be nice that the map would fit to the markers (so if I place markers only in Europe, the image would show Europe only, not the whole globe). The target display is a java webapp so I would nee to generate the image and get it in a byte form, so I can display it on the page. Does any one know an Java API that has such functionality? I looked at things like geotools or jmapviewer but they seem to lack the features I require.


texts: Meaning of ERROR 001707 when sharing layer as layer package file from ArcGIS Pro?

texts:
I have had no success when trying to share a layer as a layer package file from ArcGIS Pro 1.2.

The procedure that I use to reproduce this is:

Start ArcGIS Pro 1.2
Create a New Project using the Map.aptx template and call it Test.  I placed mine in a folder called C:\polygeo (which already exists)
In the Project pane add a folder connection to C:\polygeo where I have a point shapefile called ne_10m_populated_places.shp that is from Natural Earth
Drag and drop that shapefile into the Map
Open the Layer Properties for that new layer and use the Metadata tab to set its Description to Test Description (and then OK to close the Layer Properties)
Right-click on the layer to choose Share As Layer Package
On the Share As A Layer Package dialog:

set Save package to file (as an aside, the option to Upload package to Online account seems to work)
set Name & Location to C:\polygeo\ne_10m_populated_places.lpkx
set Summary to Test Summary
set Tags to test-tag
click Analyze (there should be no errors or warnings)
click Package

Expect to see:

  ERROR 001707: Data cannot be opened: Shape
  File(C:\polygeo\ne_10m_populated_places) ne_10m_populated_places.shp

I have not been able to find any documentation on ArcGIS Pro ERROR 001707.

Does anyone know what may be causing it?

For anyone trying to reproduce this error the ne_10m_populated_places.shp can be downloaded from http://www.naturalearthdata.com/downloads/10m-cultural-vectors/:


texts: Creating new layer in order to fix duplicate categories in GRASS?

texts:
I'm trying to run a network analysis in Grass 7.0 following this tutorial: http://www.ing.unitn.it/~grass/docs/tutorial_641_en/htdocs/esercitazione/network_analysis/node3.html

Since I'm new to Grass I already got stuck when trying to create a new layer in order to eliminate duplicate categories from my lines layer.
The query I use is:

  v.category --overwrite --verbose input=network_connect output=network_cat3 type=line layer=3 cat=1

The output is:

  ERROR: Required parameter  not set:
          (Action to be executed.)
  (Fri Aug 19 16:46:37 2016) Command executed (0 Sek)

I translated this from German, so if anything in the Error message seems fishy, that's my bad.


texts: What is the geom_type at CreateLayer() in ogr for xyz-GeoJSON?

texts:
My GeoJSON-File is a Multipolygon which creates several triangles with a X, Y and Z coordinate. Z determines the height of the specified point for X/Y.
If i try to render my GeoJSON into a shapefile by Pythan and ogr, I have to tell the CreateLayer-Function a geom_type. I've fond only Point, Line, Polygon and MultiPoint, MultiLine, and MultiPolygon.
MoltiPolygon ignores the values of height.
Is there any specification to tell the ogr in Python, that I have a MultiPolygonZ source where the height needs to be stored into the shapefile, as well?

My curent line here is:

layer = ds.CreateLayer(concenate, srs, geom_type=ogr.wkbMultiPolygon)

This works, but the resulting shapefile has no values to the height of the points.


texts: Why does Select Layer By Attribute give errors 000622 & 000628?

texts:
I try to select line segments from a layer "t_Pavement_Line which is open in TOC using the following line:

    arcpy.SelectLayerByAttribute_management("t_Pavement_Line", "NEW_SELECTION", '"ROAD_CLASS" = \'%s\'' % classname)

I get this Error message:

    ExecuteError: ERROR 000622: Failed to execute (Select Layer By Attribute). Parameters are not valid.
ERROR 000628: Cannot set input into parameter in_layer_or_view.


texts: OpenLayers3 Draw Interaction and click event issues

texts:
I plotted a cluster layer within a proximity(polygon) using draw interaction.Clusters get plotted inside.

But when i tried to put click events on cluster layer using map.forEachFeatureAtPixel , click events don't give the feature attributes of clusters within the drawn proximity even after hiding and removing the drawn polygon layer.

When i checked the layer which is clicked using map.forEachFeatureAtPixel i found it is not the cluster layer.

Please Note: Click events on clusters work fine for all cluster layers outside the drawn polygon.

Can't upload the code,sorry.


texts: Geoserver - Disallowing unauthenticated user from using getLegendGraphic Service

texts:
I created a geoserver rule as ADMIN and *.* for allowing an authenticated admin to see all services and functions.  Why does Geoserver allow Geoexplorer to retrieve the legend images for saved maps even if I haven't authenticated yet?  Is this not how service rules work?  Is this a bug? I was looking through the Geoserver JIRA and could not find any mention of this.  My current authentication chain does not allow for anonymous authentication.

I'm using Geoserver 2.5.1 as a WAR with Tomcat7.

Simple request params

http://domain.com/geoserver/wms?request=GetLegendGraphic&width=20&height=20&layer=layer&transparent=true&format=image%2Fpng&legend_options=fontAntiAliasing%3Atrue%3BfontSize%3A11%3BfontName%3AArial&SCALE=2179976.903424944


texts: Add a link to a popup window Leaflet map: adding html to setContent

texts:
I'm adding a metadata link to an info popup. I've been trying a variety of ways to do this, but all turn up errors, there appear to be different rules applying for html within setContent...

Here's the map: Click on the info button  under the zoom controls for the popup in question: https://pdxcyclesafetymap.neocities.org/

Code for the popup:

var helloPopup = L.popup().setContent("<img src='images/cog.png'style=width:150px;height:150px;
"+
"Portland Cycle Safety Map"+
"
This map is intended to bring dangerous intersections"+
"
and street segments to the attention of Portland area cyclists: this is a work in progress."+
"
Each skull marks the location a cyclist has been killed by an automobile sometime between 2005 and 2017."+
"
'metadata"'+
);


texts: Using Integrate tool in ArcGIS for Desktop on several features with range of XY Distance thresholds?

texts:
I have a feature class where there are gaps between the polygon features.  These features are supposed to be contiguous.  The most efficient way that I have discovered is using the integrate tool which allows me to select the polygon features in addition to measure the distance between the two features and use that as my XY parameter.  However I notice if I use the integrate tool on say features; a and b, then the boundary of c and b becomes incorrect based on my integrate with a parameter of X amount of feet.  My polygon layer has 23 features in it with varying XY thresh holds, how can I efficently integrate the entire polygon feature class so that my edges are contiguous for all polygons that share a neighbor?

The image below shows the "gaps".  I would select each feature where there is a gap between measure and use the measurement in the XY tolerance in the integrate tool.

I am using ArcGIS for Desktop 10 so the Aggregate Polygons tool is not available in this version.


texts: Satellite image disappearing when zooming in

texts:
I'm trying to get a satellite image using the SentinelHub plug in.

This works, however, when I zoom in to make my other layers clearer (roads etc.) the satellite image disappears and the background becomes white.

Can anyone tell me what might be causing the problem and how to fix it?

Before zooming in:

Then zooming in slightly:


texts: Choosing first letters of words in attribute table using Python Parser of ArcMap Field Calculator

texts:
I want to write a Python script in Field Calculator which can pick up or choose the fist letters of words and put them in a new field.
For example:
Lepidium catapycnon..........> Lc
Gold blossom tree............> Gbt


texts: Convert raster to smooth contour polygons in QGIS

texts:
I have had limited success creating smooth contour polygons from a raster layers.

I have created contour lines and boundary polylines and then used QGIS Union.  Then I use QGIS Polygonize.  Then use zonal statistics to the original raster to get the value on the polygons.

The 2nd method I tried was to convert raster to polygon.  Then create a new attribute column rounded to the contour I desire.  Then I dissolve lines between the same value.  Then use v.generalize.  This method is a little quicker, but it will leave some polygons jagged/unsmoothed.

Is there an easier method that actually works all the time.  Seems like there should be a simple command that does the entire process.


texts: Negative numeric values in Geoserver sql view params

texts:
I want make some query based in lon/lat values of a specific place. A simple example is this:

SELECT C.tipo, C.geom
FROM climas AS C
WHERE ST_Contains(ST_SetSRID(C.geom, 4326), ST_SetSRID(ST_Point(%x%, %y%),    4326))

As default values for x and y, I put the lon/lat coordinates of my city (-44, -3)

The problem is, when I visualize the layer without params, this works. When I put the params, if they are negative, I receive a message saying:

  error:java.lang.RuntimeException: java.io.IOException java.lang.RuntimeException: java.io.IOException java.io.IOExceptionInvalid value for parameter x

I'm using the default regular expression hinted by geoserver. Need something for negative values in the URL? Or maybe the regular expression isn't correct?


texts:
I would like to be able to specify an order of label placement priorty in QGIS. Rather than just saying 'label in the top left', I would like to be able to say 'prefer to label top left, but if not then label bottom left'.

There is already an option to be able to do this in ArcMap.

I cannot see a way to replicate this in QGIS. My current workaround is to use the offset placement from point, but with some exceptions using data defined placement coordinates. This starts to get very time consuming after a few large datasets!

Any suggestions or best practices would be welcomed here.


texts: What can I do when snapping doesn't work correctly?

texts:
I'm trying to manually edgematch some adjacent polygons. When I try to add a vertex to snap it to another vertex on an adjacent polygon, it invariably undoes itself, as shown in this video that I made demonstrating the problem (I couldn't figure out any way to post it here since it is not an image so the best I can do for now is a link).

Bad Snapping

There are 2 vertices in the beginning and I add a third. It is actually removing one of the 2 vertices on the right side and putting the result in the middle. It must be a tolerance of some sort but I tried changing it from 10 pixels to 5 (or 3 or 1) pixels and zooming way in and the problem still occurs (it just becomes more difficult to find the snap point). Is there any other setting I could modify to force the added vertices to snap without moving the other placed vertices?


texts: How can I measure(length) line of segment in postgis?

texts:
This is my first question.
I wonder that, how can I get length of line's segment?

Yellow line is my road.
blue Points are my start and end points. So I want to measure point1 and point2 between route.This points are variable. So I ll select them.
How can I get it with postgis query ?

I UPDATE MY QUESTION. :)


texts: Displaying legend item as polyline instead of polygon in ArcMap?

texts: I am using ArcMap 10.1 and am trying to create a legend that will display a polygon layer as a single line.


texts: How to automatically convert dots on paper maps to vector points?

texts:
I have several paper maps with dots that I need to convert to vector points.  I was hoping to get some suggestions on how to accomplish this without recreating each and every point by hand.  I know the basic workflow will look something like this:

Scan paper maps.
Use QGIS Georeferencer plugin to georeference maps.
Run some kind of analysis to identify dots and convert to vector
points.

I'm fine on 1-2, but a little stuck on 3.  I've provided a sample crop of one of the maps below.


texts: My program won't accept my output file path

texts:
Traceback (most recent call last):
  File "C:\Program Files (x86)\Wing IDE 101 4.1\src\debug\tserver_sandbox.py", line 21, in
  File "C:\Program Files (x86)\Wing IDE 101 4.1\src\debug\tserver_sandbox.py", line 12, in countHospitals
  File "C:\Program Files (x86)\ArcGIS\Desktop10.3\arcpy\arcpy\analysis.py", line 692, in Buffer
    raise e
arcgisscripting.ExecuteError: ERROR 000210: Cannot create output c:\Scripts\Lab 6 Datauffers.shp
Failed to execute (Buffer).

The code should perform the following tasks:

Create 3 mile buffer around the schools in the input shapefile
Count number of hospitals in each school's buffered area using Spatial Join tool
Print the name of each school with zero (0) hospitals based on the "Join_Count" field
from the resulting feature class from "Spatial Join" using search cursor

My code:

import arcpy

work = arcpy.env.workspace = raw_input("Enter the full workspace path: ")
sType = ("HIGH SCHOOL")
result = "buffer_shp"

def countHospitals(work, sType, result):

    buffering = "c:\Scripts\Lab 6 Data\buffers"
    output = "output.shp"
    # Process: Buffer
    arcpy.Buffer_analysis(output, buffers, "3 Miles", "", "", "NONE", "", "")

    # Process: Spatial Join
    arcpy.SpatialJoin_analysis(buffers, Hospitals, sJoin_shp, "JOIN_ONE_TO_ONE", "KEEP_ALL", "OBJECTID_1 \"OBJECTID_1\" true true false 10 Long 0 10 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,OBJECTID_1,-1,-1;OBJECTID_2 \"OBJECTID_2\" true true false 10 Long 0 10 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,OBJECTID_2,-1,-1;OBJECTID \"OBJECTID\" true true false 9 Long 0 9 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,OBJECTID,-1,-1;NAME \"NAME\" true true false 27 Text 0 0 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,NAME,-1,-1;ADDRESS \"ADDRESS\" true true false 28 Text 0 0 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,ADDRESS,-1,-1;FACILITY \"FACILITY\" true true false 17 Text 0 0 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,FACILITY,-1,-1;DISTRICT \"DISTRICT\" true true false 16 Text 0 0 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,DISTRICT,-1,-1;VERIFIED \"VERIFIED\" true true false 2 Text 0 0 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,VERIFIED,-1,-1;BUFF_DIST \"BUFF_DIST\" true true false 19 Double 0 0 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,BUFF_DIST,-1,-1;ORIG_FID \"ORIG_FID\" true true false 10 Long 0 10 ,First,#,C:\\Scripts\\Lab 6 Data\\buffer.shp,ORIG_FID,-1,-1;NAME_1 \"NAME_1\" true true false 33 Text 0 0 ,First,#,Hospitals,NAME,-1,-1;ADDRESS_1 \"ADDRESS_1\" true true false 28 Text 0 0 ,First,#,Hospitals,ADDRESS,-1,-1;CITY \"CITY\" true true false 7 Text 0 0 ,First,#,Hospitals,CITY,-1,-1;STATE \"STATE\" true true false 7 Text 0 0 ,First,#,Hospitals,STATE,-1,-1;ZIPCODE \"ZIPCODE\" true true false 11 Text 0 0 ,First,#,Hospitals,ZIPCODE,-1,-1;PHONE_NUMB \"PHONE_NUMB\" true true false 17 Text 0 0 ,First,#,Hospitals,PHONE_NUMB,-1,-1;AV_ADD \"AV_ADD\" true true false 28 Text 0 0 ,First,#,Hospitals,AV_ADD,-1,-1;AV_STATUS \"AV_STATUS\" true true false 1 Text 0 0 ,First,#,Hospitals,AV_STATUS,-1,-1;AV_SCORE \"AV_SCORE\" true true false 6 Long 0 6 ,First,#,Hospitals,AV_SCORE,-1,-1;AV_SIDE \"AV_SIDE\" true true false 1 Text 0 0 ,First,#,Hospitals,AV_SIDE,-1,-1;FullAddr \"FullAddr\" true true false 254 Text 0 0 ,First,#,Hospitals,FullAddr,-1,-1", "CONTAINS", "", "")

with arcpy.da.UpdateCursor(sJoin_shp, ['JOIN_COUNT', 'NAME']) as uCursor:
    for row in uCursor:
                    if row[0] == 0:
                        print ("The following schools have zero hospitals within a 3-mi radius "), row[1]
countHospitals(work, sType, result)

My code wont accept where I want to store my buffer output. I think. Any ideas as to the problem?


texts: Does HDOP correlate with GPS position accuracy?

texts:
I collected 1000+ position data from GPS receivers and realized that a position with lower HDOP does not necessarily have higher accuracy.
So my question is does HDOP correlate with GPS position accuracy?


texts: Results from split vector layer not showing

texts:
I'm trying to split a vector layer into multiple polygons (there'll be about 200 or 300, when it works), but nothing is showing up once I've run the vector split tool, neither on the layers side of the window, or in the processing results on the other side. Processing history says it's run. This is what I've done, several times:

I'm using QGIS Desktop 3.8.1 and the CSR is EPSG:4326 - WGS 84 - Geographic.

What I'm doing wrong?


texts: How to query an ArcGIS Feature Service by attribute using Python?

texts:
I have an ArcGIS Feature Server and what I have accomplished so far is this:

import json, requests

url = 'http://hebid.hessen.de/hlbk/hlbk_grundlagen/MapServer/1/query?f=json&returnGeometry=true&spatialRel=esriSpatialRelIntersects&geometry={%22xmin%22:531435,%22ymin%22:5700831,%22xmax%22:543116,%22ymax%22:5716844,%22spatialReference%22:{%22wkid%22:25832}}&geometryType=esriGeometryEnvelope&inSR=25832&outFields=*&outSR=25832'
resp = requests.get(url)
data = resp.json()
file_path = r'D:\testlayer.geojson'
with open(file_path, 'w') as f:
    json.dump(data, f)

This produces a perfect valid geojson (works with ArcMap and QGIS) using a given boundary-box, but I want to query this layer by attribute. According to this link https://resources.arcgis.com/en/help/rest/apiref/fsservicequery.html this should look something like this:

http://hebid.hessen.de/hlbk/hlbk_grundlagen/MapServer/query?layerDefs={"1":"OBJECTID=5256"}&returnGeometry=true&f=json

This does not give me my polygon as a geojson, but just the general information about the service itself. Any ideas what I am doing wrong?


texts: How to fire .NET code in my ArcGIS application from javascript?

texts:
So I'm developing a .NET desktop application for ArcGIS in VB. As a new part of the application I have developed an SVG / Javascript map as a 'loading screen' where you select your dataset to load into ArcMap. This map is currently held in a html doc and displayed in an html viewer in a windows form. My question is this:

Is it possible for javascript events/functions in my web map to fire .NET code in my application? Can I pass variables from javascript to .NET? Etc...

I imagine this will be quite tricky if possible, as you won't want any old javascript on a website to start messing with your running .NET applications. I'm willing to hack this all together if there's a way.

Thanks.

EDIT: To be clear, this is not an ASP.NET (server side) web application, this is a purely vb.net desktop application. So I think AJAX / Postback solutions won't work in this case.


texts: About speed, commitChanges too slow. How avoid this?

texts:
I am working on a standalone script, executed out of QGIS console,  and I have noticed that the script takes a lot of time when it executes the line: layerB.commitChanges ()

How my script works ?

I have two files that are identical, one is the copy of another.

The script loads both files (fileA and fileB), then it searches for similar points (very near, distance <= 15m. ) over fileA,  and deletes it over the fileB. When finished, do commit over the fileB

What would be the fastest way to do this task?
 I hava read this-> Is it possible to change the geometry without startEditing the layer
But I do not finish understanding the proposed solution.

The input data table is approximately 200,000 points and the output table is almost half.

layerA = QgsVectorLayer("fileA.sqlite", "layerA", "ogr")
layerB = QgsVectorLayer("fileB.sqlite", "layerB", "ogr")

featA = [feat for feat in layerA.getFeatures()]
featB = [feat for feat in layerB.getFeatures()]

del_feat = []

target_spatial_index = QgsSpatialIndex()
# populate the spatial index
for f in featB:
    target_spatial_index.insertFeature(f)

with edit(layerB):
    for feat in featA:

        # Skip if into deleted list
        if feat.id() not in del_feat:

            point1 = feat.geometry().asPoint()

            nearest_ids = target_spatial_index.nearestNeighbor(point1,2)

            high_id=0
            for id in nearest_ids:

                outFeat = layerB.getFeatures(QgsFeatureRequest().setFilterFid(id)).next()

                #skip himself
                if outFeat.id() != feat.id():

                    # distance between points
                    point2 = outFeat.geometry().asPoint()
                    dist = sqrt(point1.sqrDist(point2))

                    high_id=outFeat.id()

                    if dist <= 15:
                        # remove from spatial index
                        target_spatial_index.deleteFeature(outFeat)

                        layerB.deleteFeatures([high_id])
                        # add to list of deleted
                        del_feat.append(high_id)

layerB.commitChanges()
crs=QgsCoordinateReferenceSystem("epsg:32629")
print "Saving file..."
# Save memory layer to file
error = QgsVectorFileWriter.writeAsVectorFormat(layerB, "fileB.sqlite", "UTF-8", crs , "spatialite")


texts: How to prepare domain on multifeature using Arcmap?

texts: I am using a Personal Geodatabase.  How do I set up domains on multiple fields?


texts: Possible for client to edit geoJSON on server?

texts:
I am a webmapping newbie so pardon any misnomers I am probably making.

I have a geoJSON file that represents garbage cans as lines running along a street (i.e. the line is an aggregated representation of all the cans on that street, so if 1 can is 20 lbs and there are 3 cans on the street, the line feature on that street will have a weight attribute of 60 lbs).

These lines have four attributes: weight, time (seconds per picking up all cans on that street), number of cans on the street, and the route number that services them. These attributes appear in a div at the bottom of the screen when you click/select a given line:

Is it possible to allow the user/client to change the route attribute value in a manner that it changes the original geoJSON and then refresh the web page to see the new changes (lines are symbolized by route number)?


texts:
I tried to use the same coding from Loading Multiple CSV Files into QGIS 2.10 to load multiples CSV in QGIS but in QGIS 3.6 and it does not work. The code provided in the above link was used and the following error shows:

  Traceback (most recent call last):
    File "C:\PROGRA~1\QGIS3~1.6\apps\Python37\lib\code.py", line 90, in runcode
      exec(code, self.locals)
    File "", line 1, in
    File "", line 13, in
  NameError: name 'QgsMapLayerRegistry' is not defined

I am new to python so I don't really know what is going on. =P

Maybe some changes in parameters or so? Was some function deprecated from Python 2.7 to 3?

How the answer from here would need to look like now? To facilitate, this is the code:

import glob, os

# Define path to  directory of your csv files
path_to_csv = "C:/Users/You/Desktop/csv folder/"

# Set current directory to path of csv files
os.chdir(path_to_csv)

# Find each .csv file and load them as vector layers
for fname in glob.glob("*.csv"):
uri = "file:///" + path_to_csv + fname + "? delimiter=%s&crs=epsg:4326&xField=%s&yField=%s" % (",", "x", "y")
name = fname.replace('.csv', '')
lyr = QgsVectorLayer(uri, name, 'delimitedtext')
QgsMapLayerRegistry.instance().addMapLayer(lyr)


texts: Is there a way to project LiDAR data from a GCS to a state plane coordinate system?

texts:
I have LiDAR data that is in GCS NAD83. I have converted that LiDAR data into an LAS Dataset for use in ArcPro. I'm now trying to convert that LAS dataset into a raster and the tool keeps failing with the the message: Error 000210: Cannot create output . Results window below.

ESRI's documentation states that this error results from either a data lock, incorrect path or limited access rights. I know this isn't the case, however. The output directory is on the C drive of my pc. I have also tried exporting to a tif file and gotten the same results. I've come across others that have had this issue and it seems that the problem frequently occurs if the original LiDAR data is in a GCS. Using lastools, I was able to confirm that our data is in NAD83 GCS. How can I convert that to a projected coordinate system in order to successfully run this tool?


texts: Measuring distance in spherical Mercator vs zoned UTM

texts:
I've got points in WGS84 lat/long and I'd like to measure "small" (less than say 5km) distances between them.

I can use the haversine formula from http://www.movable-type.co.uk/scripts/latlong.html and it works very well.

I'd like to use Python Shapely libraries though, so that I can do more operations than just distance, and because at the scale I'm working with, a flat earth is a good enough approximation.  To reliably project the geographic coords to a cartesian coord, I'm using Python's proj4, but seem to get bigger errors than I'd like.

If I use the local UTM zone, I get differences between haversine of a couple of meters, which is fine.  But I don't want to have to work out the UTM zone (the points could be worldwide), so I tried with "spherical Mercator" but now the differences between haversine and projected distances are well over 100%.  Is this really right for spherical Mercator?  All I really want is a workable Cartesian projection for two points within 5km of each other anywhere in the world.

from shapely.geometry import Point
from pyproj import Proj

proj = Proj(proj='utm',zone=27,ellps='WGS84')
#proj = Proj(init="epsg:3785")  # spherical mercator, should work anywhere...

point1_geo = (-21.9309694, 64.1455718)
point2_geo = (-21.9372481, 64.1478206)
point1 = proj(point1_geo[0], point1_geo[1])
point2 = proj(point2_geo[0], point2_geo[1])

point1_cart = Point(point1)
point2_cart = Point(point2)

print "p1-p2 (haversine)", hdistance(point1_geo, point2_geo)
print "p1-p2 (cartesian)", point1_cart.distance(point2_cart)

At this point, the haversine distance between them is 394m, and using utm zone 27, 395m.  But if I use spherical Mercator, the Cartesian distance is 904m, which is way off.


texts: Changing the label position of one item in a layer at QGIS

texts:
I am working on a mapping project of Australia and need to display the layer labels respectively to other project items.
So, I have placed the labels in the center quadrant of the polygon.

And then I pushed them done the Y axis by three millimeters

However, I want to push done only the label Of Western Australia and no other label from this layer.
Is it passable to change the position of only one label from the same layer?


texts: Looking for a good online map with rivers

texts:
Is there any internet map source (zoomable, tiled, i.e. "google like") with good display of rivers? I.e. well labelled even on larger scales? So that I'm for example able to track Donau river (with label) on scale like this.

Basic google map is just terrible, rivers are not well displayed, often lost on larger scales, and in 90% cases not labelled.
OpenStreetMap - similar, furthermore they denote roads as blue (!!!) which is even more confusing! (really bad practice)
OpenCycleMap - even worse...


texts: Selecting rectangle box from raster tiled table in PostGIS

texts:
I am new to PostGIS as well as SQL queries for rasters.

I have imported a raster .tif file into PostGIS 2.1 with tiles. I would need to select a bounded rectangle from within this raster using SQL statements but i have no clue on how to begin.

Can i ask what are the possible/better ways to do that?

Do i need to have the top left and bottom right coordinates/lat long to get the rectangle?

I had search and read many other related posts about using ST_Intersects, BBOX and ST_MakeEnvelope but i don't really quite get them.


texts: Mapbox.js: Accessing data markers created in Mapbox Studio

texts:
I used Mapbox Studio to create some projects because my data set was too large to use Mapbox Editor. I loaded my data into a source project and then used that source as a layer in a style project, which I then uploaded to Mapbox here.

I used conditional formatting rules in CartoCSS to color-code the data markers based on the "species" represented (the map is created from migration data of large marine animals). I would like to create another map using Mapbox.js that allows me to change the marker color based on whether that species is selected or not, like this one.

If a species is not selected, the marker should be dull gray in color, if the species is selected, it will be its vibrant color.

My question is: how do I reference those individual species markers that have been cooked as a source project into a layer of the style project in Mapbox Studio? Can I use Mapbox.js to access the CartoCSS or will I need to access some GeoJSON? The data was originally uploaded as csv so I would need to understand the GeoJSON working behind the Mapbox scenes (if that is what is happening).


texts: CORS ERROR with own WMS GeoServer layers

texts:
I am building a weather application. For this I want to use OpenLayers as the base for my application and then I want to be able to add my own made layers from PostGIS with GeoServer to this application.

At the moment I have GeoServer running and made a map with OpenLayers using Node.js, NPM and Parcel.

When I go to localhost:1234 (where this application is shown) I can see the OpenLayers basemap, but I won't show the WMS-layer from GeoServer. I get the following message:

Access to image at 'http://localhost:8080/geoserver/....' from origin 'http://localhost:1234' has been blocked by CORS policy: No 'Access-Control-Allow-Origin' header is present on the requested resource.

I have searched all over for the solution, but I think I am missing something. I have already installed CORS on NPM, added CORS to the GeoServer web.xml and the Tomcat web.xml. (Below the code I added.)

Still I get this error. I have been searching for days to solve this, but I just cannot find what I am missing and it is starting to get really frustrating. Anybody ideas?

EDIT:
I found that in de code for my application there was the following line of code:

crossOrigin: 'anonymous'

I removed this code, because I am not working with crossOrigin but with CorsFilter in my other code. Now the console displays:

Cross-Origin Read Blocking (CORB) blocked cross-origin response http://localhost:8080/geoserver/.... with MIME type text/xml. See https://www.chromestatus.com/feature/5629709824032768 for more details.

END EDIT///

Code added in Tomcat web.xml:

>   <filter>
>     <filter-name>CorsFilter</filter-name>
>     <filter-class>org.apache.catalina.filters.CorsFilter</filter-class>
>   </filter>   <filter-mapping>        <filter-name>CorsFilter</filter-name>
>       <url-pattern>/*</url-pattern>   </filter-mapping>
>     <init-param>
>         <param-name>cors.allowed.origins</param-name>
>         <param-value>http://localhost</param-value>
>     </init-param>
>     <init-param>
>         <param-name>cors.allowed.methods</param-name>
>         <param-value>GET,POST,PUT,DELETE,HEAD,OPTIONS</param-value>
>     </init-param>
>     <init-param>
>         <param-name>cors.allowed.headers</param-name>
>         <param-value>Content-Type,X-Requested-With,Accept,Authorization,Origin,Access-Control-Request-Method,Access-Control-Request-Headers</param-value>
>     </init-param>
>     <init-param>
>         <param-name>cors.exposed.headers</param-name>
>         <param-value>Access-Control-Allow-Origin,Access-Control-Allow-Credentials</param-value>
>     </init-param>

Code added to GeoServer web.xml

<filter>
    <filter-name>CorsFilter</filter-name>
    <filter-class>org.apache.catalina.filters.CorsFilter</filter-class>
    </filter>
    <filter-mapping>
        <filter-name>CorsFilter</filter-name>
        <url-pattern>/*</url-pattern>
    </filter-mapping>


texts:
I'm having trouble getting the ListDatasets (or any of the list functions) to return anything but an empty list when pointed at our SDE geodatabase, but the same code works for a file geodatabase. I searched how to connect to the SDE GDB and I can't see a mistake. I know the connection file works, I can open it with Catalog and log in just fine and it has a saved password with it. I'm at a loss, which means I have probably done something stupid. Help?

# Import arcpy module
import arcpy
from arcpy import env

env.workspace = r"F:Replicated Data\Auto_Extract_Connection.sde"

datasetList = arcpy.ListDatasets()

print datasetList

env.workspace =  r"F:\Replicated Data\Gatekeeper_Extract.gdb"

datasetList = arcpy.ListDatasets()

print datasetList

Returns:

[]

[u'ElectricDataset', u'CPUDataset', u'TransmissionDataset']


texts:
I am working on a project to create a centerline from a polygon automatically. The user select a polygon and click a button, the centerline will create for him.

I have tried the build-in function “Collapse Dual Lines to Centerline”, but it cannot perform well so I decide not use this function in this project. In further study, I find out some terms such as “Medial Axis”, “Skeleton”, “Voronoi polygon”, those are some concepts related to centerline formation. However, I find difficult to understand those algorithm.

Does anyone have a suggestion to me for this project?
Does anyone have a suggestion for how to implement those concepts in ArcMap?
Is there any sample script or algorithm I can refer to as reference?

Development language: C#, VB.Net
ArcGIS : ArcMap 9.3.1


texts: Removing vector layer for overwriting using PyQGIS

texts:
In QGIS, I've tried to remove a layer, execute a processing and overwrite the input layer. This in my Python code:
from qgis.core import *

for lyr in QgsProject.instance().mapLayers().values():
    if lyr.name() == "test":
        QgsProject.instance().removeMapLayers([lyr.id()])

infn =  QgsProject.instance().homePath()+"/temp.shp"
outfn = QgsProject.instance().homePath()+"/temp.shp"

processing.run("native:refactorfields", {'INPUT':infn,\
            'FIELDS_MAPPING':[{'expression': '\"Field1\"','length': 10,'name': 'Field1','precision': 0,'type': 10},\
            {'expression': '\"Field2\"','length': 50,'name': 'Field2','precision': 0,'type': 10}],\
            'OUTPUT':outfn})

The layer is successfully removed from the legend but I received this error for the output file:
error OGR: C:/temp/temp.shp is not a directory.

It would seem as if the layer remained in memory and did not allow to overwrite.
Any help?


texts: How to create an SQL Persisted Computed Column which values are derived from spatial column

texts:
Our GDB is using the SQL native GEOMETRY data storage type (SQL 2012 Standard) and the spatial reference for this uses a unit of square feet. I want to create a (persisted) computed column which will populate another column as acres.

Using the following query I'm able to produce the results I want:

Select shape.STArea()/43560 as Area from DATAOWNER.CTF_STANDS

In setting up the computed column I used the following:

Editing on the front end GIS application (ArcGIS), the Shape.STArea() field updates during edits, but updates to my computed column fail to happen.

I think that Shape.STArea is actually a function and appears to be enforced/maintained by the SDE/GDB Schema itself, so working with spatial column appears to be the best way forward.

I don't want to convert the Geometry column itself to another unit, but rather use it to populate another column via a formula which converts square feet to acres.

How do I create a computed column which will derive acres from the geometry column who's current values are in square feet?

This other SE post's answer is very similar.


texts: Converting a geocentric latitude to a geodetic latitude

texts:
I'm working with a 3D model of the earth (WGS84).

The input to my problem is a geocentric lat-long. Is there an efficient way to convert this to a geodetic lat-long?

My initial attempt to solve this seems to work correctly: First, convert the geocentric lat-long back to it's corresponding point in 3D space. Then, feed that XYZ into a library function that returns the geodetic latitude.

It's straightforward, but the conversion from XYZ to geodetic lat-long is not cheap because it uses some fixed number of iterations to come up with a reasonably accurate result... I am looking for something that runs faster.

Given the geocentric latitude, is there a more efficient route I can take to get the corresponding geodetic latitude?


texts: Stacking lidar-derived height raster as ancillary data for segmentation in ENVI?

texts:
I get a a black and white image after following these steps.  Does anyone know why I'm not getting a color aerial image beneath the mask?

My setup consists of these steps:

I mosaicked several NAIP tiles together (ArcGIS Pro: Mosiac to New Raster), then clipped that raster to a polygon representing my study area.
I created a DTM (ground) and DSM (first return) from a LAS dataset (ArcGIS Pro: LAS Dataset to Raster) using the study area polygon as the processing extent, NAD83 UTM 13N as output coord sys. Units are in meters, and I applied a Z factor of 3.2808399 to convert to feet.
I used the Minus tool to substract the DTM from the DSM to get a height band.
I created an NDVI polygon to use as a mask in ENVI.
I used the Segment Only Feature Extraction tool. Input: NAIP raster. Mask: NDVI polygon. Ancillary data: height band. Custom bands: NDVI.

In ENVI, here is the NAIP raster clipped to the study area:

Here is the height band raster:

Here is the ndvi mask:

Here is the stacking result from step 5:

There should be an aerial image beneath the mask, not a black and white image. I tried stacking without using the height band as ancillary data, and it worked the way it should. I have used this method successfully to stack bands for segmentation before. It seems obvious there is something wrong with the height band, but I cannot figure out what it is!

Has anyone encountered this before?

Using ArcGIS Pro 2.3.2, ENVI 5.5 v3.2.


texts: Automating ArcGIS append Geoprocessing Tool

texts: I have about 100 feature classes inside twenty feature datasets inside a geodatabase. We have about 50 individual geodatabases. I want to load all geodatabases one by one in a one single master geodatabase. I believe anyone with experience in Python programing can help me. I don't have any coding experience.


texts: Image quality and aliasing effects in GeoWebCache WMTS tiles in QGIS

texts:
I am accessing WMTS data in QGIS 3.8 (and QGIS 3.4) from a GeoServer at the following address:

https://gis.lmi.is/geoserver/gwc/service/wmts?REQUEST=GetCapabilities

The layer is LMI_Kort and the WMTS tiles are served via a GeoWebCache.  The data appear distorted/aliased, with zig-zag effects on straight lines e.g. roads.  It happens in the native projection and if reprojected.  I have tried setting the zoom level (via Scale box on GUI) to match the levels for which the tiles were created e.g. 1:10,000, 1:25,000, 1:50,000 etc., but it doesn't help.

There are no issues when the same layer is viewed directly from the WMS server, however WMS is significantly slower:

https://gis.lmi.is/geoserver/ows?version=1.3.0

There is some discussion of a similar issue on the QGIS mailing list (http://osgeo-org.1560.x6.nabble.com/WMS-GeoWebCache-image-quality-in-QGIS-td5157450.html) but it is unresolved.

All the WMS/WMTS options are set to their defaults.  Is there a way to configure QGIS to render the WMTS tiles without distortion?

Update:

I tried using the TMS URL to add the layer as XYZ tiles:

https://gis.lmi.is/geoserver/gwc/service/tms/1.0.0/LMI_Kort@EPSG%3A3857@png/{z}/{x}/{y}.png

This brings in tiles in strange locations.  If I use {-y}, to invert the y axis, I don't get any tiles at all.

Update:

Using different DPI-Mode settings when adding the WMTS connection in QGIS doesn't help. (all, off, QGIS, UNM, GeoServer)

Update:

Using the different scales provided in the demonstration maps on the GeoServer e.g. 1:273k doesn't help.

Update:

The effect is the same on both a hi-dpi laptop screen and a normal desktop monitor.


texts: Replacing each point in a point feature with a specific shaped polygon

texts:
I have many random points in one layer. I want to replace each individual point with a specific shaped polygon like the image below illustrates.

Kind of like applying some sort of irregular buffer to each point.

The reason for this is that I want to identify where symbols that have the shape of the polygon I intend to use overlap.

Is this possible? I would like to apply this to several thousand points at once also. Maybe with python, model builder or FME?


texts: Can't generate true orthophoto from my software

texts:
I am a GIS specialist and a surveyor. Lately I'm trying to use as a surveying tool my DJI Phantom 3 pro combining with different software for generating DEM, DSM, 3d model and orthophoto. But I can't generate true orthophoto from the software I used like Agisoft or pix4d, orthophoto after generated is deformed like in the pictures below.

My workflow was:

I added the photos
I Loaded the camera positions (with EXIF info on)
I aligned the photos
I optimized camera alignment with these settings

Camera accuracy: 10
Marker accuracy: 0.005 (not utilized)
Scale bar accuracy: 0.001 (not utilized)
Projection accuracy: 0.1
Tie point accuracy: 4

I builded Dense Cloud with quality Medium and depth filtering Aggressive
I builded mesh with surface type Height field, source data Dense cloud and face count Medium, Interpolation enabled
I builded texture with mapping mode Orthophoto, blending mode Mosaic, texture size 8192
Then I exported as orthophoto with blending mode Mosaic, enable color correction off

How can I generate a true orthophoto from this, what I did wrong here?


texts: How to import projection from a file in gdalwarp?

texts:
I want to batch clipping and transforming rasters by a shp file. It looks like the projection file importing is wrong:

import subprocess
import os
import gdal,osr

inshape = r'D:\NorthGrass\Map\NorthTemperateGrass2.shp'
inpath = r'D:\GIMMS'
outpath = r'D:\NorthGrass\GIMMS82_132'

inproj = r'D:\NorthGrass\Map\GCS_WGS_1984.prj'
outproj = r'D:\NorthGrass\Map\WGS_1984_Albers.prj'

sr_proj = osr.SpatialReference()
sr_proj.ImportFromWkt(inproj)
tr_proj = osr.SpatialReference()
tr_proj.ImportFromWkt(outproj)

for filename in os.listdir(inpath):
     if os.path.splitext(filename)[1] == '.tif':

          inraster = os.path.join(inpath,filename)
          print filename+'   begin processing'

          outraster=os.path.join(outpath,filename)

          subprocess.call(['gdalwarp','-s_srs',sr_proj,'-t_srs',tr_proj,'-cutline',inshape, '-of', 'GTiff',\
                 '-crop_to_cutline','-tr','8000','8000','-ot','Int16',\
                 '-dstnodata','-32768',inraster, outraster])
          print outraster+'   is done'
print ' ----All Done----'

Hint error:

TypeError: argument of type 'SpatialReference' is not iterable


texts:
How do I make geoserver update layers and stores? I've updated the tables in PostGIS used by a data store and layer, but I cannot seem to get it to update. I have the new data running, but I had to create a new Store and Layer just to get it to work. Is there an easier way to refresh stores and/or layers?

Thanks


texts: Reversing polyline direction based on raster value using ArcPy?

texts:
I am trying to create an ArcPy script that will read a raster cell value at the start point and end point of a polyline, and flip the line based on those values.

I do not want to create another feature class, I just want to edit a feature class in an existing file geodatabase.

I am able to get the cell values at the nodes, however, the lines are not flipping direction, so I think the problem lies within the second for loop.

Is there a way to do this without writing the geometries to a list and reversing the order because I have had more trouble trying to do just that?

The code is as follows:

import arcpy

inFeatures = "Lines"
shapeName = arcpy.Describe(inFeatures).shapeFieldName
rows = arcpy.SearchCursor(inFeatures)

Raster = 'srtm2f.tif'

for row in rows:
    feat = row.getValue(shapeName)
    xy1 = feat.firstPoint
    xy2 = feat.lastPoint
    resultSTART = arcpy.GetCellValue_management(Raster, str((xy1)), "1")
    CellValueSTART = float(resultSTART.getOutput(0))
    resultEND = arcpy.GetCellValue_management(Raster, str((xy2)), "1")
    CellValueEND = float(resultEND.getOutput(0))

for features in inFeatures:
    if CellValueSTART < CellValueEND:
         arcpy.FlipLine_edit(inFeatures)
    elif CellValueSTART > CellValueEND:
        print 'Feature was correctly digitized'
    elif CellValueSTART == CellValueEND:
        print 'Feature had same elevation'


texts: How to connect GeoServer with Java to a PostGIS/PostGres database?

texts:
I am having trouble using GeoServer to perform a shortest path analysis between two network points. My question is about connecting GeoServer with Java to a PostGIS/PostGres database.

I know that in MapServer we can use PHP functions like PG_CONNECT to connect and PG_QUERY to execute SQL queries... but what are the equivalent commands for GeoServer?


texts:
I have 2 feature class files of contour lines for an area of interest - from 2 different sources. We would like to select just one to use and I'm looking for any methods to effectively compare the files.

I'm using ArcGIS as my primary application to run an analysis.

Wondering if converting the lines to surfaces and comparing the surfaces is appropriate, or comparing values only along the lines itself?


texts: wine + LAStools in QGIS 2.8.1, Ubuntu 14.04

texts:
I can't get wine + LAStools to work in QGIS 2.8.1 (ubuntu 14.04). I only get the 22 open source functions that comes with QGIS. Outside qgis wine + LAStools works fine.

My guess is that I am not specifying the folders correctly under Processing Options -> Tools for LiDAR data, and I can't find any documentation about it. You can enter the LAStools folder and the wine folder but it is not clear how this should be done. I looked at the source code for LAStoolsUtils.py where I believe this is partly handled. I think this line is important.

folder = wine_folder + "/wine " + lastools_folder

My LAStools folder can be found at:

/home/.../.wine/drive_c/LAStools/QGIS_toolbox/processing

So I put:

[wine folder] /home/.../

[LAStools folder] /drive_c/LAStools/QGIS_toolbox/processing

I thought /.wine could be the problem but changing /wine to /.wine did not help. I also tried all kinds of variation on the directory paths but nothing worked.

Any ideas how to get this to work? There are other questions like this but for windows/MAC. E.g.
Using LAStools in QGIS gives "global name 'os' is not defined"?


texts: I wonder if there is a possibility (plugin etc.) to create a summarising diagram of attributes more or less comparable to the summarising statistics. e.g. having a shape with a high number of polygons of different landuse attributes and want to create a diagram with the sum of area for the different land use types as a bar chart or pie chart diagram. x axes land use types, y axes sum of area.


texts:
I created a script in which I run numerous algorithms from processing. Each resulting layer is stored in a shp file. I just want to keep the final layer and delete the others, but when I want to erase them while QGIS is still open, they are still in usewhile not shown in the layers list, I can't erase them. This is a problem when I run my script a 2nd time because it doesn't work.

I tried removing the layers with
QgsMapLayerRegistry.instance().addMapLayer(layer)
QgsMapLayerRegistry.instance().removeMapLayer(layer.id())
but I had no success so far...
.prj .qpj and .shx are deleted, only remains .dbf and .shp.

Here is the code

import processing
from qgis.core import *

layer=processing.getObject("carte")
layer.invertSelection()

processing.runalg('qgis:fixeddistancebuffer', layer, "20","5",False, "c:/Users/LENOVO/Desktop/zones/tampon+.shp")
processing.runalg('qgis:dissolve', "c:/Users/LENOVO/Desktop/zones/tampon+.shp", True, None, "c:/Users/LENOVO/Desktop/zones/dissolve.shp")

tampon=QgsVectorLayer("c:/Users/LENOVO/Desktop/zones/tampon+.shp", 'finale','ogr')
QgsMapLayerRegistry.instance().addMapLayer(tampon)
QgsMapLayerRegistry.instance().removeMapLayer(tampon.id())

processing.runalg('qgis:fixeddistancebuffer', "c:/Users/LENOVO/Desktop/zones/dissolve.shp", "-20","5",False, "c:/Users/LENOVO/Desktop/zones/carte_finale.shp")
dissolve=QgsVectorLayer("c:/Users/LENOVO/Desktop/zones/dissolve.shp", 'finale','ogr')

QgsMapLayerRegistry.instance().addMapLayer(dissolve)
QgsMapLayerRegistry.instance().removeMapLayer(dissolve.id())
layer = QgsVectorLayer("c:/Users/LENOVO/Desktop/zones/carte_finale.shp", 'finale','ogr')
QgsMapLayerRegistry.instance().addMapLayer(layer)


texts: visualising edge_graph from list pgr_dijkstra graph

texts:
I am creating a query that finds the shortest path from dangerous crossing nodes to hospitals (many-to-many). I have the proper output from pgr_dijkstra, with output as follows

seq | path_seq | start_vid | end_vid | node | edge | cost | agg_cost
-----+----------+-----------+---------+------+------+------+----------
   1 |        1 |         2 |       3 |    2 |    4 |    1 |        0
   2 |        2 |         2 |       3 |    5 |    8 |    1 |        1
   3 |        3 |         2 |       3 |    6 |    9 |    1 |        2
   4 |        4 |         2 |       3 |    9 |   16 |    1 |        3
   5 |        5 |         2 |       3 |    4 |    3 |    1 |        4
   6 |        6 |         2 |       3 |    3 |   -1 |    0 |        5
   7 |        1 |         2 |       5 |    2 |    4 |    1 |        0
   8 |        2 |         2 |       5 |    5 |   -1 |    0 |        1

I am trying to use st_makeline to create seperate lines for each individual node path, however thus far the function has been creating one big line
SELECT  st_transform(st_makeline(t1.the_geom),4326)
        FROM pgr_dijkstra(
                'SELECT osm_id as id, source, target, cost FROM spatial.melbourne_edges',
                array(select hospital_id from hospitals_in_vicinity),
                array(select crossing_id from dangerous_crossing_nodes), false) as route
                inner join spatial.melbourne_nodes t1 on route.node = t1.id --);

How am i able to make a column of all the different paths?


texts: Create elevation dependent mean values in ArcGIS RasterCalculator

texts:
I have two rasters. Raster1 is a reclassified raster for certain elevation value ranges. Raster2 is an elevation raster too but with gaps. I want to fill the gaps of Raster2 by the median of the elevation class of Raster1. But Raster1 is only used to give the class areas. So I want to calculate the median of Raster1 taking the elevation classes/areas of Raster2 and save them to the gaps of Raster1.

The theoretical idea is more or less clear but I don't know how to create the raster calculator statement.


texts:
Okay, my goal is to mask a geotiff by means of a shapefile.

Firstly, I made sure that the shapefile and geotiff had same projections:

from osgeo import ogr, osr
import gdal
import os

driver = ogr.GetDriverByName('ESRI Shapefile')
inDataSet = driver.Open(r'shapefile.shp')

# from Layer
inLayer = inDataSet.GetLayer()
spatialRef = inLayer.GetSpatialRef()
# input SpatialReference from Geometry
inFeature = inLayer.GetNextFeature()
geom = inFeature.GetGeometryRef()
inSpatialRef = geom.GetSpatialReference()

raster_path = 'raster.tif'

# output SpatialReference
raster = gdal.Open(raster_path, gdal.GA_ReadOnly)
wktproj = raster.GetProjectionRef()
outSpatialRef = osr.SpatialReference()
outSpatialRef.ImportFromWkt(wktproj)

# create the CoordinateTransformation
coordTrans = osr.CoordinateTransformation(inSpatialRef, outSpatialRef)

# create the output layer
outputShapefile = r'shapefile_out.sh'
if os.path.exists(outputShapefile):
    driver.DeleteDataSource(outputShapefile)
outDataSet = driver.CreateDataSource(outputShapefile)
outLayer = outDataSet.CreateLayer("GSHHS_f_L1_out_layer", geom_type=ogr.wkbMultiPolygon)

# add fields
inLayerDefn = inLayer.GetLayerDefn()
for i in range(0, inLayerDefn.GetFieldCount()):
    fieldDefn = inLayerDefn.GetFieldDefn(i)
    outLayer.CreateField(fieldDefn)

# get the output layer's feature definition
outLayerDefn = outLayer.GetLayerDefn()

# loop through the input features
while inFeature:
    # get the input geometry
    geom = inFeature.GetGeometryRef()
    # reproject the geometry
    geom.Transform(coordTrans)
    # create a new feature
    outFeature = ogr.Feature(outLayerDefn)
    # set the geometry and attribute
    outFeature.SetGeometry(geom)
    for i in range(0, outLayerDefn.GetFieldCount()):
        outFeature.SetField(outLayerDefn.GetFieldDefn(i).GetNameRef(), inFeature.GetField(i))
    # add the feature to the shapefile
    outLayer.CreateFeature(outFeature)
    # dereference the features and get the next input feature
    outFeature = None
    inFeature = inLayer.GetNextFeature()

# Save and close the shapefiles
inDataSet = None
outDataSet = None

Okay, now I use the outputtet shapefile which should have the same projection as raster now to mask:

import fiona
import rasterio
import rasterio.mask
import gdal

img_name = 'raster.tif'
shape_name = r'shapefile_out.shp'

raster = gdal.Open(img_name, gdal.GA_ReadOnly)
wktproj = raster.GetProjectionRef()

with fiona.open(shape_name, "r") as shapefile:
    geoms = [feature["geometry"] for feature in shapefile]

with rasterio.open(img_name, driver='GTiff', crs = 'wktproj') as src:
    out_image, out_transform = rasterio.mask.mask(src, geoms, crop=False, invert=True)
    out_meta = src.meta.copy()

out_meta.update({"driver": "GTiff",
                 "height": out_image.shape[1],
                 "width": out_image.shape[2],
                 "transform": out_transform})

with rasterio.open("masked.tif", "w", **out_meta) as dest:
    dest.write(out_image)

But I get this result??

What's going on?


texts: Finding obsolete/dated scripts from retired arcscripts.esri.com?

texts:
How can I find scripts/tools/add-ins which are linked from any website and reference the retired http://arcscripts.esri.com/ site?

I first try looking them up in the new site, https://codesharing.arcgis.com but most of the time I cannot find what I am looking for in https://codesharing.arcgis.com and try to find cached (or archived or other links referenced) versions of the source to at least find the publisher, so I can do an in-depth web search.


texts: I'm looking at the possibility of using a kml layer in a 'geoserver' to provide attribute information on a point. Is there a KML file for the Wikipedia Places articles as it is used in Google earth or it is a Google 'locked' functionality?


texts: QGIS difference operation with XYZ tiles

texts: I have a regular layer from a database containing WKT polygons. I'm also using a XYZ Tiles layer of Google Maps. I would like to perform the Vector > Geoprocessing Tools > Difference operation between the two layers, except using an XYZ Tiles layer is not a possibility in the Difference dialog box (notice no expand arrow to allow access to OpenStreetMap or gmm (Google Maps)). QGIS 3.16 being used.


texts:
I am looking for a good interactive mapping software/suite that will allow me to display polygons over a topography map, and be able to display HTML when a polygon is clicked, for more information. I have looked into Google Maps, but the problem is that this map will be in a static location, without internet.

Goal:  Find  a way to make interactive maps offline.

Researched: Google Offline caching, but from what i understand, it will not be able to cache an entire country with polygons, and HTML files. I have found GMapCatcher, can anyone verify it caches very well, and in full detail?

Other than that, any suggestions are welcome.

This map will be an interactive map for a WW II exhibit, which will consist of polygons( describing battles and positions), and when clicked on, they are to display some information about the object.


texts: leaflet: Different information popups for each polygon

texts:
I have a geoJSON file of all the counties in Michigan. I have embedded this file in a Leaflet map.

I want each of these county polygons to display the name of the county in a pop-up when clicked on. This is basically a duplicate of Different popup for every polygon (leaflet), except that my map has 83 polygons, and I want to find a way to do this without calling L.bindPopup() 83 times in a row.

Should I be looking into some kind of javascript function or for loop?

I am not very familiar with javascript at all.

My javascript code thus far:

var mymap = L.map('mapid').setView([43.293218, -84.607686], 7);
L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png', { attribution: ' OpenStreetMap' }).addTo(mymap);

function addDataToMap(data, map) {
    var dataLayer = L.geoJson(data);
    dataLayer.addTo(map);
}

$.getJSON("county_coords.geojson", function(data)
    {
        addDataToMap(data, mymap);
    });

Here is a sample of my geoJSON file (county_coords.geojson). Maybe there is a way I can call the "LABEL" property in the geoJSON file to display the county name in a popup?

{
  "type": "FeatureCollection",
  "crs": {
    "type": "name",
    "properties": {
      "name": "urn:ogc:def:crs:OGC:1.3:CRS84"
    }
  },
  "features": [
    {
      "type": "Feature",
      "properties": {
        "FIPSCODE": "001",
        "FIPSNUM": 1,
        "NAME": "Alcona",
        "LABEL": "Alcona County",
        "TYPE": "County",
        "CNTY_CODE": "001",
        "SQKM": 1798.557217,
        "SQMILES": 694.423636,
        "ACRES": 444428.3814,
        "VER": "13a",
        "LAYOUT": "landscape",
        "PENINSULA": "lower",
        "URL": "http://per510a.mtri.org:8090/download/attachments/327710/NC_001_Alcona_2005.zip?version=2&modificationDate=1182876202960&api=v2"
      },
      "geometry": {
        "type": "MultiPolygon",
        "coordinates": [
          [
            [
              [
                -83.886988854009,
                44.540800273944
              ]...