texts: Use self-defined Leaflet identifiers to work with jsTree checkboxes

texts:
I have some hierachic structured point-layers in a GeoJSON file using jsTree (ajax) with checkboxes to add them to the map. As I have about 20.000 points structured in 2.000 nodes representing the layers I have to to it in an automated way.

My idea is to connect the jsTree id with an id in the GeoJSON, so the Node ID of the jstree corresponds with the "treeid" (node.id===treeid, both are strings) in the GeoJSON file. The string ID's are generated in a Database and exported as JSON. Each checkbox adds an ID assigned layer to a layergroup on the map. Unchecking should remove that particular layer. The user can activate as much checkboxes as he likes.

Based on this question I set the default _leaflet_id equal to the "treeid" to use it as an identifier for each Layer. This works basically.

Example of the GeoJSON:

var collection = {"type": "FeatureCollection","crs": {"type": "name","properties": {"name": "urn:ogc:def:crs:OGC:1.3:CRS84"}},   "features": [{ "type": "Feature",
   "properties": {
     "dta_1": "Node_1",
     "treeid": '01 01 002 00001'
   },
   "geometry": {
     "type": "Point",
     "coordinates": [38.933333, 40.816667]
   }
 },
 {
   "type": "Feature",
   "properties": {
     "dta_1": "Node_2",
     "treeid": '01 01 002 00002'
   },
   "geometry": {
     "type": "Point",
     "coordinates": [38.9551, 40.7626]
   }}]}

Example of the jsTree (simple JSON):

$('#data').jstree({'core': {'data': [{
     "id": "01",
     "parent": "#",
     "text": "Root"
   }, {
     "id": "01 01 002 00001",
     "parent": "01",
     "text": "Node_1"
   }, {
     "id": "01 01 002 00002",
     "parent": "01",
     "text": "Node_2"
   } ...

My checkbox code:

var group = {};
var layergroup =  L.layerGroup().addTo(map);
$('#data').on("check_node.jstree", function(e, data) {
            group =  L.geoJson(collection, {
            onEachFeature: function(feature, layer) {
                           layer._leaflet_id = data.node.id;},
            filter: function(feature, layer) {
                    if (feature.properties.treeid==data.node.id) return true}});
            group.treeid = data.node.id
            group.addTo(layergroup);
 /*         layergroup.eachLayer(function(layer){
            layer._leaflet_id = data.node.id;
            }) */
     })
.on("uncheck_node.jstree", function(e, data){
layergroup.eachLayer(function(e){
if (e._layers[data.node.id]) layergroup.removeLayer(e);})});
//   if (layer._leaflet_id === data.node.id) {
//   layergroup.removeLayer(layer)}

A working fiddle is here: jsfiddle

Using the "onEachFeature" Function works fine for one point, it is removed correctly. So I thought, it works, but as you can see in the fiddle, if there are many points (Node_2) again only one point (the first one) is added to the map. (Because only the first one is assigned?)

Using the "eachLayer" method adds all points as expected, but of course the uncheck code doesn't work. Changing to the outcommented code does not work.

Where is the point I'm missing on the "onEachfeature" Function to add all points with the same treeid?

Or the other option: How to change the uncheck code?

I need to add all points with the same id referencing a specific node.id to a layergroup by a checkbox, unchecking it should remove them from the layergroup and the map.


texts: Creating Concentric circles using the ggmap function from GGMAP package R

texts:
I would like to create two concentric circles using the ggmap package and ggmap function, one being 1/4 mile and the other 1 mile radius around the 16 metro rail stations in Los Angeles California. It looks like ggmap doesn't scale the size in geom_point(aes(size=1)) to miles, so I would like to know how to convert this to miles. (I looked at the dismo package but I would like to use ggmap to do this.)

library(ggmap)

d = matrix(c(-118.2342,  34.05620,
       -118.2461,   34.05490,
       -118.2513,   34.04932,
       -118.2587,   34.04863,
       -118.2749,   34.05637,
       -118.2901,   34.06270,
       -118.2917,   34.07653,
       -118.2917,   34.08991,
       -118.2918,   34.09771,
       -118.3081,   34.10174,
       -118.3252,   34.10163,
       -118.3385,   34.10155,
       -118.3627,   34.14000,
       -118.3768,   34.16850,
       -118.3015,   34.06175,
       -118.3088,   34.06169), 16, 2, byrow=T)

colnames(d) = c("long", "lat")
d = data.frame(d)

m = get_map(location = "Hollywood,California", maptype = "roadmap", zoom = 12)
la_map = ggmap(m)
# HERE I WANT TO ADD A RADIUS OF 1 MILE and A RADIUS OF 1/4 MILE OVER THAT.
la_map +
geom_point(aes(x = long, y = lat, size = 1, col = "blue"), data = d) +
geom_point(aes(x = long, y = lat, size = 1/4 ,col = "red"), data = d)


texts: Layer missing from Export to PDF from ArcPy?

texts:
I am using the arcpy module to loop through a number of map features (community areas), zoom to them, title the map with the name of the community area, and export to pdf. The trouble is that while the title, legend, and line features export, a symbolized polygon layer doesn't. It still shows up in the legend, but the map doesn't show it.

Trying a couple of sanity checks, I have

confirmed that that layer in question is visible in the ArcMap map that I'm starting with,
tried making sure that that layer is on top,
instructed my code to make that layer the only visible one before exporting, and
had Python list the set of layers that are visible, right before exporting.

None of those tests give me anything better, or give me any insight to the problem.

My code is below. Since I'm also a relative newbie to coding with arcpy, any suggestions that you have for solutions or more elegant coding would be very welcome!

My code is:

import arcpy, os
myPath    = r"C:\Users\<myname>\Documents\Heat Maps\\"
myMapPath = r"C:\Users\<myname>\Documents\Heat Maps\Map Output\\"

arcpy.env.workspace = myPath
myMap = arcpy.mapping.MapDocument(myPath + "Heat Maps of Need - v4-0 - Added Streets for Zoomed in Maps.mxd")

AllLayers = arcpy.mapping.ListLayers(myMap)
df = arcpy.mapping.ListDataFrames(myMap)[0]
lyrList = ["Needs Index"]
nComm = 77

# Clear visibility of all layers except for the community areas layer

for lyr in AllLayers:
    if lyr.name = "Needs Index":
        print "Turning on Layer " + lyr.name
        lyr.visible = True
    else:
        lyr.visible = False

# Loop through specified layers and export to pdf

if os.path.exists(myPath + "Heat Maps of Need Indicators.pdf"):
    os.remove(myPath + "Heat Maps of Need Indicators.pdf")
PDFdoc = arcpy.mapping.PDFDocumentCreate(myPath + "Heat Maps of Need Indicators.pdf")

for l in lyrList:

    # Turn on proper heat layer
    for lyr in AllLayers:
        if lyr.name == l:
            lyr.visible = True
    arcpy.RefreshActiveView()

    # Check which layers are visible
    for lyr in AllLayers:
        if True == lyr.visible:
            print "Layer " + lyr.name + " is visible"

    # Loop through all community areas
    for c in range(1, nComm + 1): # Add 1, since the top value of the specified range is not included

        # Run tasks specific to the community area in question
        for lyr in AllLayers:
            if lyr.name == "Community Areas":

                # Get name of community area
                SqlWhere = " \"AREA_NUMBE\" = \'" + str(c) + "\'"

                CommunityRowObj = arcpy.SearchCursor(lyr, SqlWhere)
                for row in arcpy.SearchCursor(lyr, SqlWhere):
                    CommunityName = row.getValue("COMMUNITY")
                    CommunityName = CommunityName.title()

                # Zoom in to community area
                arcpy.SelectLayerByAttribute_management(lyr, "NEW_SELECTION", SqlWhere)
                df.zoomToSelectedFeatures()
                arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")

        for x in arcpy.mapping.ListLayoutElements(myMap):
            if x.name == "Title":
                x.text = "Heat Map of Needs Indicators\n" + l + " - " + CommunityName

        arcpy.RefreshActiveView()

        # Export the layer to PDF
        arcpy.mapping.ExportToPDF(myMap, myMapPath + "Needs Heat Map - " + l + " - " + CommunityName + ".pdf")
        PDFdoc.appendPages(myMapPath + "Needs Heat Map - " + l  + " - " + CommunityName + ".pdf")

    # Turn off current Needs Layer to clear it from the next one
    for lyr in AllLayers:
        if lyr.name == l:
            lyr.visible = False

arcpy.RefreshActiveView()
arcpy.RefreshTOC()
PDFdoc.saveAndClose()

del PDFdoc
del AllLayers


texts: Semi automatic classification plugin, downloading sentinel 2, authentication problem

texts:
When I try to find a Sentinel 2 image for download with this plug in, it asks me about my credentials, I enter them, but nothing happens.

On the ESA page, I can login, no problem.

I use OSX 10.8.5 and QGIS 2.14.8

How do you download Sentinel 2 imagery with this plug in? Why doesn't it work for me?
Can I download the images directly from the ESA page and then upload them into QGIS and the plug in?


texts: Combining Bathymetric and Elevation Data using R?

texts:
Here is some elevation data of some islands.

I wish to add depth info to the water. With the generous help of @Spacedman at Overlaying elevation data and bathymetry data in R? I get:

To resolve the blockyness of the lower resolution water data I tried, resample(to_harbour, to_island, method = "bilinear") yielding:

I'm not sure if that is correct. Note that I have lost some of the water around the islands. I'm not sure why this is nor am I sure how I should identify the water in to_island. Currently I am just looking at the data and guessing that values less than 75.109m are water and should be looked up from the bathymetric raster layer.

Dropbox:

to_island

to_harbour

Code

library(sp)
library(raster)

library(tidyverse)
library(rayshader)

to_island <- readRDS(file = "to_island.RDS") # https://www.dropbox.com/s/hc86gzdt6nm4j4w/to_island.RDS?dl=0
to_harbour <- readRDS(file = "to_harbour.RDS") # https://www.dropbox.com/s/awl4clyb3254n14/to_harbour.RDS?dl=0

# Create matrix for 3D plotting
to_island_matrix = matrix(raster::extract(to_island, raster::extent(to_island), buffer = 10000),
                          nrow = ncol(to_island), ncol = nrow(to_island))

# 3D plot land without water depth info ####################################################
to_island_matrix %>%
  sphere_shade(zscale = 0.1, texture = "imhof1") %>%  # Lower zscale maps points to wider range of colours
  plot_3d(to_island_matrix, zscale = 1, fov = 0, theta = 170, phi = 30, windowsize = c(1000, 800), zoom = 0.35,
          water = TRUE, waterdepth = 75.2, wateralpha = 0.5, watercolor = "#88DDFF",
          waterlinecolor = "white", waterlinealpha = 0.5)

render_snapshot()
rgl::clear3d()

###########################################################################################
# Add water depth values to land raster layer ##############################################

# Resample to_harbour to match higher resolution of to_island
to_harbour_res <- resample(to_harbour, to_island, method = "bilinear")

# Identify points on to_island that are water
w <- which(to_island[] < 75.109)
length(w)

# Get coordinates of those points (water)
xy <- SpatialPoints(to_island)[w]

# Add CRS
proj4string(xy) = projection(to_island)

# These locations aren't in the coordinate system of the water raster so we need to transform them (WGS84 v NAD83)
xyt = spTransform(xy, projection(to_harbour))

# Get the water raster values at those points:
ex = raster::extract(to_harbour_res, xyt)

# Replace the locations in the land raster with the values sampled from from the water raster:
to_island[w] = ex

# Create matrix for 3D plotting
to_island_matrix = matrix(raster::extract(to_island, raster::extent(to_island), buffer = 10000),
                          nrow = ncol(to_island), ncol = nrow(to_island))

# # 3D plot land WITH water depth info
to_island_matrix %>%
  sphere_shade(zscale = 0.1, texture = "imhof1") %>%  # Lower zscale maps points to wider range of colours
  plot_3d(to_island_matrix, zscale = 1, fov = 0, theta = 170, phi = 30, windowsize = c(1000, 800), zoom = 0.35,
          water = TRUE, waterdepth = 75.2, wateralpha = 0.5, watercolor = "#88DDFF",
          waterlinecolor = "white", waterlinealpha = 0.5)

render_snapshot()
rgl::clear3d()


texts:
Given a FeatureCollection that contains multiple Features (each consisting of polygons), I would like to use the Python API to loop through them all, and batch export them to Google Cloud Storage.

However, I can't figure out how to make all the code server-side, and hence my code runs extremely slowly (using the major no-no of getInfo inside a loop).

for i in range(0, 1000):
  region = ee.Feature(myFeatCol.get(i)).geometry()

  task = ee.batch.Export.image.toCloudStorage(
    image = naip.clip(region),
    description = 'uniqueName'+str(i),
    fileNamePrefix = 'uniqueName'+str(i),
    bucket = 'myGoogleBucketName',
    scale = 1,
    region = region.getInfo()['coordinates'])

  task.start()

The objective here is to end up with 1000 images in the bucket. It works, but is painstakingly slow! Several-minutes-per-image slow.

How can I create an efficient server-side loop to export an entire FeatureCollection?


texts: Print Multiple PDFs Using Different Layers Using a QGIS 3.0 Layout

texts:
I have an open QGIS 3.0 project with a layout called "plan_view1".  There are numerous selected layers (>100) that I want to turn on one-by-one and print the "plan_view1" layout to a pdf.
This is my code:
def printpdfmulti(layoutname):
    selected_layers = qgis.utils.iface.layerTreeView().selectedLayers()
    projectInstance = QgsProject.instance()
    layoutmanager = projectInstance.layoutManager()
    layout = layoutmanager.layoutByName(layoutname) #Layout nameprojectInstance = QgsProject.instance()
    for layer in selected_layers:
        qgis.utils.iface.legendInterface().setLayerVisible(layer, True)
        exporter = QgsLayoutExporter(layout)
        exporter.exportToPdf("C://data//" + layer.name() + ".pdf", QgsLayoutExporter.PdfExportSettings() )
        qgis.utils.iface.legendInterface().setLayerVisible(layer, False)

printpdfmulti("plane_view1")

The error messsage I get is:
Traceback (most recent call last):
  File "C:\OSGEO4~1\apps\Python36\lib\code.py", line 91, in runcode
    exec(code, self.locals)
  File "<input>", line 1, in <module>
  File "<string>", line 13, in <module>
  File "<string>", line 7, in printpdfmulti
AttributeError: 'QgisInterface' object has no attribute 'legendInterface'

What have I done wrong?
I know my python programming is still very rough on the edges.


texts: GDALOpen: "ERROR 1: Too Big for Int" is informative only?

texts:
When I use GDALOpen to open an ESRI Imaging .img file I see the error:
"ERROR 1: Too Big for Int"
Is this an actual error, or a warning? I see no way to avoid the message, It is puzzling since I am only opening the file, not trying to read from it. I seem to be able to read the file properly anyway. Can anyone explain this?


texts: ArcGIS Server - Change Virtual Directory

texts:
I recently installed ArcGIS for Server 10.3.

Is there a way for me to change the virtual directory for my map services?

For example, instead of using

http://mywebsite.com/arcgis/rest/services/

I want to use

http://mywebsite.com/103/


texts: Calculating mean anual precipitation on GEE

texts:
I am using GEE to calculate mean anual precipitation between 2000 and 2019 using
GPM: Global Precipitation Measurement (GPM) V6. Because there are 20 years in the series, I first summarize all 20 years in a new image collection "ee.ImageCollection(precip_mensual)" which contains mean precipitation for all 20 Januaries, all 20 Februaries, and so on (steps 1 & 2).
In the resulting collection, each pixel has 12 values (one for each month). Then, next I need to somehow summarize that into an average year (step 3).

Using the "Inspector", I am not getting the values I am supposed to in points for which I know the mean value for anual precipitation. That being said, I think that either think my approach is not the correct one, or it is but there is something missing in my script.

1) Add shapefile of study area M (using asset)

Map.addLayer(M,{},'M');

var GPM = ee.ImageCollection('NASA/GPM_L3/IMERG_V06')
          .filter(ee.Filter.date("2000-06-01", "2019-06-30"))
          .select('precipitationCal');

2) Inter-anual aggregation
Here, I use a function to derive the average year at a monthly scale.

var months = ee.List.sequence(1, 12);
var precip_mensual = months.map(function(m) {
var precip_men = GPM.filter(ee.Filter.calendarRange(m, m, 'month')).mean();
return precip_men;
});

var precip_mensual = ee.ImageCollection(precip_mensual);
print('Precipitacion', precip_mensual)

3) Intra-anual aggregation

var precip_mean = precip_mensual.mean().rename('precip-med').clip(M)
Map.addLayer(precip_mean, {}, 'Precipitation-Mean')

Map.centerObject(M, 4)

4) Create region

var ExportArea
=ee.Geometry.Rectangle([-77.9881000000000029, -55.9185000000000016,
-37.4998999999999967, 2.0680000000000001]);
Map.addLayer(ExportArea, {color: 'FF0000'}, 'poly');

Export.image.toDrive({
image: precip_mean,
description:'Precipitation_mean',
folder:'WildBoar',
crs:'EPSG:4326',
maxPixels: 1e13,
scale:1000,
region:ExportArea,});


texts: Javascript API: How can I change renderer on the fly with feature layer hosted on ArcGIS Online?

texts:
I am new to Javascript and somewhat new to GIS (especially web applications of GIS). I'm trying to understand how I can have a drop down menu or html buttons to change which data is rendered in my map.

I have a feature layer hosted on ArcGIS Online, so I don't think I can use ArcGISDynamicMapServiceLayer. This is only if you have ArcGIS Server, right? I only have ArcGIS Desktop.

I can, however, customize my renderer using ClassBreaksRenderer. I thought it would make sense to set up a few different renderers, like so:

            var approvalRate = new ClassBreaksRenderer(symbol, "MEAN_appro");
                approvalRate.addBreak(0,.19999, new SimpleFillSymbol().setColor(new Color([255, 0, 0, 0.5])));
                approvalRate.addBreak(.2,.39999, new SimpleFillSymbol().setColor(new Color([255, 128, 0, 0.5])));
                approvalRate.addBreak(.4,.59999, new SimpleFillSymbol().setColor(new Color([255, 255, 0, 0.5])));
                approvalRate.addBreak(.6,.79999, new SimpleFillSymbol().setColor(new Color([139, 209, 0, 0.5])));
                approvalRate.addBreak(.8,1, new SimpleFillSymbol().setColor(new Color([56, 168, 0, 0.5])));

            var averageMortgageAmount = new ClassBreaksRenderer(symbol, "MEAN_dolla");
                averageMortgageAmount.addBreak(0, 24.9999, new SimpleFillSymbol().setColor(new Color([255, 0, 0, 0.5])));
                averageMortgageAmount.addBreak(25,49.9999, new SimpleFillSymbol().setColor(new Color([255, 128, 0, 0.5])));
                averageMortgageAmount.addBreak(50,99.9999, new SimpleFillSymbol().setColor(new Color([255, 255, 0, 0.5])));
                averageMortgageAmount.addBreak(100,199.9999, new SimpleFillSymbol().setColor(new Color([139, 209, 0, 0.5])));
                averageMortgageAmount.addBreak(200,Infinity, new SimpleFillSymbol().setColor(new Color([56, 168, 0, 0.5])));

So, I have two different renderers set up, called approvalRate and averageMortgageAmount.

I thought I could set up some logic to change which renderer is applied based on user input. So I have set up the renderer variable, get its value from a button, and use that to determine which renderer to use with a switch. Like below.

            var renderer;

            function getRenderer(rendererFromButton) {
                renderer = rendererFromButton;
            }

            switch(renderer) {
                case "Approval Rate":
                    featureLayer.setRenderer(approvalRate);
                    break;
                case "Average Mortgage Amount":
                    featureLayer.setRenderer(averageMortgageAmount);
                    break;
                default:
                    featureLayer.setRenderer(approvalRate);
            }

This ends up with a map that correctly renders approvalRate, but never actually changes dynamically. The buttons work right (I've tested the output of the variable "renderer" value after pressing the buttons and it is correct).

It seems like the problem lies in the browser calling the map server again. I don't know what I'm talking about here, but I assume once the map is rendered by the server and sent to the browser as a graphic layer or whatever, it can't be called again and re-rendered.

Is that true?

Here's what I have so far: http://matthewiannowlin.com/planning/thesis/

This is long-winded, but I figured more detail gets a more helpful answer.

Bottom line: is there anyway I can change which field is rendered from an ArcGIS Online feature layer with user input? Do I need to link to different pages that each have different renderers defined? Or do I need to just get ArcGIS Server?


texts: How to open existing sqlite DB in QspatiaLite plugin?

texts:
Using QspatiaLite plugin in QGIS I'm able to create a new database and work with it (load/edit/import/export data). But how can I add a connection to an existing database?
How to add an existing database to the drop-down-menu of the available database connections?


texts:
I have a shapefile of lines containing an amount n of elements. I exported these features and analyzed them with another program (Matlab), which provided me as an output a row vector of length n containing only 1s and 0s. I would like to read this vector in ArcMap, and highlight somehow the elements whose Objtect ID is equal to a vector elements equal to 1.

The reason why I want to highlight them is to export an image showing the chosen features.

I still do not have a code since I do not know if what I want to do can actually be done, nevertheless, here there is a non-python, Matlab-like code that hopefully will better describe what I would like to do:

load Portfolio.mat %this is my n-vector of 1s and 0s
for i=1:length(Portfolio)
    if Portfolio(i) ==1
        highlight (feature(ObjectID == i)) %this is the function that I do not know if it exists
    end
end

Can this be done? Is there a geoprocessing tool able to do this or is it necessary to write a python script for this purpose?


texts: Converting set of points to solid object?

texts:
I have a data generated as the output of some numerical processes. I need them to be converted into 3D solid objects (as CAD solid format) to make them importable for some particular numerical analysis software such as FEM. One initial stage is to do some Boolean operations on them. If they were is solid format it is OK, but the data is only coordinates and that's all!

The data set consist of points in which a set of points makes a particular shape say a cube. Each object has own associated points. If we draw lines between points of each object a proper object appears. These data need to be exported in 3D format consisting of only solid objects. Then a series of Boolean operations will apply.

So the questions are:

How to convert such a data set to for example a dxf format consisting only solid objects?

In particular purpose, what is the structure (format or so) of a solid object in CAD so we may write some simple code to convert our data set to solids?

(Note that the data set is large so we need finally to do some programming in any language required.)


texts: Adding break/vertex on a polyline from points

texts:
I have a polyline which consists of just one geometry. I have points which are evenly spaced (25m). How would I add these point vertices to the line so that the points will create one line segment/polyline per point vertex?

This is so that I can have an attribute for the polyline each 25m so I can use it for labeling.

I am using QGIS.


texts:
As I long ago figured out that I need more complex query-asking abilities (other than the standard select provided by arcgis), and as I will someday have to admin a small scale SDE database that will import data from my corporate's big and country-wide SDE, I wanted to set up a small workgroup SDE on a standalone computer so I could deal with it a bit, to get the hang of dealing with a sde.

But where do I start? I can't find a good manual to start with - The Arcgis 9.3.1 Help(on which I'm working) isn't much help. I know I need to Install SQL server, but I can't understand any of the other steps (creation, import of data and such).

Do you know of any good guides and resources to guide me with the arcgis 9.3.1 SDE from scratch?

(Oh, and does SDE 10 differ a lot from SDE 9.3.1?)


texts: Split Line at intersection and attach attributes

texts:
I want to slit line at intersection and attach attributes and export
I have roads table with columns id,road_name,road_type,geom with 2115 records. I try using this query - select  st_astext((st_dump(st_union(geom))).geom) from roads;
it successfully splits the line gives 5114 records but when i am trying to attach attributes it fails select  id,road_name,road_type,st_astext((st_dump(st_union(geom))).geom) from roads;

I need to split line at intersection also keep attributes means line 1 divide in 3 having same id is 1 , road_name is abc similar to line 2 divide in 2 having same id is 2 , road_name is xyz and so on.


texts: QGIS3 Shape Digitizing Circles is Drawing Improperly

texts: I'm trying to draw a circle in QGIS 3 with the shape digitizing toolbar. I've tried several of the "add circle" options and they all look normal while drawing, but once I click to finish the shape it is no longer a circle. Please see the attached photo. Any ideas on how to fix this so I can create circles?


texts: PosGIS KNN returns different results from QGIS hub spoke function

texts:
Following example here I have tried to get nearest neighbours where code is;

CREATE TABLE model.hub_spoke_spokes3f AS
SELECT
   hub.id AS hub_id,
   spoke.id,
   spoke.point_geom,
   ST_Distance(hub.geom, spoke.point_geom) AS dist
FROM hub_data.geom_data_unique AS hub
JOIN LATERAL (
SELECT id,
     point_geom
FROM spoke_data.data_set
WHERE ST_DWithin(hub.geom, point_geom, 150000) --geom is 277000, units metres
AND hub.id <> id
ORDER BY hub.geom <-> point_geom
LIMIT 100
) AS spoke
ON TRUE;

I then compared the results with QGIS Hub and spoke function which I was trying to replicate.

As screen shot below the results differ. The lines are created by QGIS but the spokes (smaller dots) are from the query above. Whats causing the difference?


texts: QGIS Server doesn't render any fonts on Ubuntu 16.04

texts:
I just updated qgis-server package to latest 3.2.x from Ubuntugis repository (ubuntugis-unstable) and everything works fine - but font rendering. All labels on tiles, legend and print templates show up as squares:

After upgrading I got a lot of this messages in nginx error log:

QFontDatabase: Cannot find font directory /usr/lib/x86_64-linux-gnu/fonts - is Qt installed correctly?

I could fix that by adding a symbolic link from /usr/lib/x86_64-linux-gnu/fonts to /usr/share/fonts and updating font cache with fc-cache -fv. But this brought no changes to visualization of labels.

The QGIS project uses Ubuntu font installed from ttf-ubuntu-font-family package and shows up fine on QGIS Desktop. The font is also listed on the server using fc-list.

There are no error messages related to fonts neither in nginx error log nor in QGIS Server logging. Could be related to Qt so, any hints?


texts: Geographic coordinates system to Lambert Conformal Conic projected CRS using a transformation matrix

texts:
Given a LCC projection with two standard parallels, is it possible to define a transformation matrix to map any given lat/lon in the covered area into the projected CRS (e.g. similar to a rotation matrix or other affine transformations).
Since the coordinates transformation equation coefficients (e.g.  theta, rho) depend on the input coordinates, I guess it's not possible to exactly calculate such a matrix, but I was wondering if any numerical approximation exists.


texts: Restoring original Settings Menu in QGIS customization?

texts:
Making some tests in Windows XP for customizing a user interface in QGIS 1.8 (see http://www.duif.net/qgis/doc/output/html/en/docs/user_manual/introduction/qgis_configuration.html#customization), I have unchecked, among other options, the Settings Menu. This means that, after restarting QGIS, there's no way of recovering the older settings because there is not already a Setting Menu. So, I cannot reload any ini file.

For your information, as a precaution, previously I have saved in a folder both the original ini file, with all the options, and the customized ini file with only some options (of course, not included here the necessary Settings Menu).

Before reinstalling QGIS (version 1,8), I would like to try other smarter options if there's any. Where and what does QGIS look at where this applications starts?


texts: Changing locale on a per app basis in ArcGIS Online

texts:
Our organisation in ArcGIS Online is localised to a non-english language. I'm now working on a project aimed at english speaking visitors.
As such, I want to change the localisation for the Web Map Application for this specific project.

Is it possible to change the localisation of the AGOL GUI for a specific web map app, without changing the localisation for our entire organisation?


texts: how to permanently add raster to satellite image in qgis

texts:
As  mentioned in the title i need to integrate a raster to a satellite image but have no idea as to how to do this.
I tried the add grid function and it works nice but is not permanent. I insert it via map composer but as this inhibits zooming I need to save it on the image. The person who will use the image as a map wishes to be able to zoom just as in any image.
How is it possibly to permanently add the a raster with 10 square kilometre grid size permanently to a satellite image I have as TIFF format?


texts: Reading single row in PyShp with shapefile.Reader

texts:
I use PyShp and shapefile.Reader to read shapefile files.
If I want to read a whole file I use:
with shapefile.Reader(file) as shp:
    for row_num, row in enumerate(shp.iterRecords()):
        print(row)

But I would like to read a single row from shapefile file with shapefile.Reader, for example row number 100. I read documentation and I didn't find any solution. How can I do that?


texts: Why do I get a shattered-glass effect when I turn on hillshade for my IDW surface?

texts:
I made an interpolated IDW surface with a cell size of 500. In the symbology tab of ArcGIS, there is the "Use hillshade effect" checkbox underneath the classification area. When unchecked, my surface looks like this:

After I check the "Use hillshade effect" box, it looks like this:

With the hillshade effect turned on, why does it look like there is a shattered-glass looking mess underneath the surface? It doesn't seem to correspond to the class breaks. What can I do to get rid of that without messing with cell size? I set the IDW surface's cell size to 500 so that the surface would look smoother/less pixelated, but the lower the cell size, the more pronounced the shattered-glass effect looks.


texts: How to I limit the x and y axes with tmap in r

texts:
I have made this plot with the following code:

library(tmap)
tm_shape(RH2) +
      tm_fill(col="AREAL") +
      tm_borders()

But there is a dot on the far right (an island) which I don't understand why it is there. Is there a way to limit the plot to avoid all the white spaces? like coord_fixed(xlim=c(7, 12), ylim=c(52, 58)) in ggmap?


texts:
I have about 300 records in an attribute table that each need an image to go with them. I've created a raster field on the table, and I've got all the needed images as rasters in the same GDB. I've been going through each record and manually adding the appropriate raster in editor.

Is there a script out there for doing this all at once?


texts:
I am trying to print a map like the image below. And I have a problem with grid and graticule.
I understand that the best way to have a grid (in UTM, meters) and a graticule (in GMS) is by adding both grids and changin the SRC to each one.

The problem with this method is that I can't find an option to add the padding marked in the photo.

The other way that I was trying is adding two Maps to the composer, the first one with the layers, and with the UTM grid in blue, and a second Map canvas backwards, in blank, showing the GMS graticule in WGS84.

I don't know if there is a better option to compose a simmilar map.

Does anybody know which would be the best method?


texts: Disable wheel scroll zooming with ArcGIS Javascript API 4.11

texts:
I want to disable wheel scroll zooming so the viewer can scroll down the page without losing the cursor in the map.

But this code freezes scrolling altogether, so when the cursor is over the map, there is no scrolling at all.

view.on("mouse-wheel", function(event) {
          // disable mouse wheel scroll zooming on the view
          event.stopPropagation();
        });

Codepen is https://codepen.io/nadinetrahan/full/VNYxRg


texts:
I am wondering if there is a way to apply a translation in coordinates for a GeoJSON export.  Currently, there are negative coordinate values in my export, which I cannot have.  I am trying to use this GeoJSON data to create a map in javascript, and I consulted the developers of the javascript library and was told that if my x values were translated 160, making them positive, the map would display correctly.  Is there an easy way to translate these values in the export or possibly do this with an expression?  I am very new to QGIS so I am unsure how to approach this problem.

Basically, I want to change values like this

[ -160.07380334546815, 22.004177347957729 ],

to this

[ 0.07380334546815, 22.004177347957729 ],


texts: Cannot find cloned python arcgispro-py3 environment with conda env list

texts:
I am trying to clone my python environment that comes with ArcGIS Pro (arcgispro-py3 located in C:\Program Files\ArcGIS\Pro\bin\Python\envs\arcgispro-py3) to a location I can access (I do not have administrative privileges). I followed these instructions from ESRI https://support.esri.com/en/technical-article/000020560 but gave a direct path name for the clone (C:\Users\me\AppData\Local\ESRI\conda\envs\arcgispro-py3_clone).
When I cloned arcgispro-py3, I did so within my Python Command prompt. The cloned environment shows up when I list out my environments in the Python Command prompt with conda env list. However, when I do so from the Anaconda prompt, the environment is not shown (only base shows). I can, however, activate it using the path to the cloned environment. (Though oddly when I run conda env list again from the cloned environment, it has no output.)
I had used this same workflow on a different computer previously and had no problems. My main question is why can I not see all of the environments with conda env list? The only thing that I can think of is that I also installed Anaconda (even though conda comes with ArcGIS Pro). I don't know if that confused anything.


texts: Adding button to QGIS toolbar by using Python and Macros in Project Properties

texts:
I've successfully created a script that adds a button to the plugin toolbar in QGIS, by running the script in the python console.

Here is my code:

import webbrowser
icon_path = MyPath

def open_url():
webbrowser.open(MyURL)

action = QAction(name)
action.triggered.connect(open_url)
action.setIcon(QIcon(icon_path))
iface.addToolBarIcon(action)

I have read that you can make QGIS to add the button when a certain QGIS project is loaded by using 'Project properties' -> 'Macros'.

Though, This won't work for me.

What I have tried is...:

import webbrowser
from qgis.PyQt.QtGui import *
from qgis.PyQt.QtWidgets import *
from qgis.utils import iface

def openProject():

url = 'https://google.dk'

def show_url():
    iface.messageBar().pushMessage()
    webbrowser.open(url, new=1)

action = QAction('Show Time')
action.triggered.connect(show_url)
action.setIcon(QIcon(MyPath)
iface.addToolBarIcon(action)

Can anybode figure out what I am missing in the code for make the Macro work out well?


texts: Error message when running sDNA license manager, Windows 10

texts:
I want to install sDNA for QGIS (Windows 10, QGIS 3) and when I try to run the license manager, I get this error message (I'm running both as admin)

I've tried to click continue, but after I copy and paste my license in the field and click on "Get Challenge URL", the same message pops up again. I clicked on "continue" again, but the status information field just says "Please wait" and nothing happens.  Reinstall didn't help and using different versions of sDNA (3.4.6 and 3.4.7) didn't help either (both times the same problem)

Do you maybe know what went wrong and what I can do to solve this problem?


texts: How can I style polygons with parameters inside the FeatureCollection properties on Leaflet?

texts:
I have a GeoJSON file that has both the FeatureCollection and the Feature types inside it. Basically my code is the following (working sample here):
L.geoJSON(myJSON, { style: style }).addTo(mymap);

function style(featureObject) {
    return {
        fillColor: getColor(featureObject.properties.name),
        weight: 2,
        opacity: 1,
        color: 'white',
        dashArray: '3',
        fillOpacity: 0.7
    };
}

function getColor(name) {
    return name === "name1" ? "#CF8562" :
            name === "name2" ? "#FFEDA0" :
                                "#000000";
}

And a small sample of the structure of my GeoJSON is the following:
[
  {
    "type": "FeatureCollection",
    "properties": {
        "name" : "name1"
    },
    "features": [
      {
        "type": "Feature",
        "properties": {
            "randomParameter" : "text"
        },
        "geometry": {
          "type": "Polygon",
          "coordinates": [
            [
            [
            -46.6479847244904,
            -23.553060554172923
            ],
            [
            -46.6479847244904,
            -23.553060554172923
            ],
            [
            -46.64805562776844,
            -23.55318890961171
            ],
            [
            -46.64826795788915,
            -23.552928264599316
            ],
            [
            -46.6479847244904,
            -23.553060554172923
            ]
            ]
          ]
        }
      }]
  },
  {
    "type": "FeatureCollection",
    "properties": {
        "name" : "name2"
    },
    "features": [
      {
        "type": "Feature",
        "properties": {
            "randomParameter" : "text2"
        },
        "geometry": {
          "type": "Polygon",
          "coordinates": [
            [
            [
            -46.648050596629034,
            -23.55393832474017
            ],
            [
            -46.64758222900779,
            -23.554141824100373
            ],
            [
            -46.64767437978837,
            -23.554322319586415
            ],
            [
            -46.64814729501603,
            -23.55411425749883
            ],
            [
            -46.648050596629034,
            -23.55393832474017
            ]
            ]
          ]
        }
      }]
  }
]

This code is not going to work properly because it doesn't recognize the featureObject.properties.name value. This value exists only inside the FeatureCollection and not the Feature type in my GeoJSON. So, I'm not capable of coloring the polygons according to the data that I have inside the GeoJSON.
I'd like to edit the style of all elements under a FeatureCollection type. However, when I print the object arriving on the style function, I see that it only sees the elements with the Feature type. I've read answers like this one telling that we can edit the style of the FeatureCollection type but I didn't manage to make it work.
I'd like to solve this problem without having to change my GeoJSON structure or changing the properties inside FeatureCollection to Feature. Is it possible? Can I make the style function change the style of everything under the FeatureCollection type instead of the Feature type in my GeoJSON?


texts: How to combine 2 vector shapes created from raster files

texts:
In Qgis 2.12 Lyon, using windows 10

I downloaded 2 raster files from an area of my interest. The rasters were on information about precipitation and evapotranspiration. Each cell had a value. I wanted to combine the information from both. It means: Calculate for each cell, the difference between precipitation and evapotranspiration.
I didn't know why, but couldn't find the way to do it in raster files. So I thought, let's convert these files in vector polygon files, so each cell now becomes a polygon with a value. My intention was to have the attribute table with all the cells and do the difference with the calculator.

But the thing is, when I polygonized (converted the raster into vectors) Some of the cells, are grouped, so 3 cells form a polygon, making it difficul for me to create a combination, or a calculation, cell-by-cell.

See images;

These were the images for precipitation info. Look what happens with the ET:

My wish would be to combine both vector shapes, and calculate the difference for each cell, but if the polygons are each one different.... How could I do it?


texts: Obscure import GeoJSON to PostGIS by ogr2ogr

texts:
I import UTF-8 encoded GeoJSON files and get a different column value in PostgreSQL.

Execution command one template is:

 ogr2ogr -f "PostgreSQL" PG:"host=127.0.0.1 dbname=db user=user
 port=5432 password=1" "file.geojson" -nln kk1

Examples:

GeoJSON file:

 ..."features": [ { "type": "Feature", "properties": { "ID": 4420590,
 "NAME": "", "CAD_NUM": "05:19:000044" }, "geometry": { "type":
 "Polygon", "coordinates": [ [ [ ...

In PostgreSQL column (cad_num) table value is:

  05:19:44

Other GeoJSON file

 ..."features": [ { "type": "Feature", "properties": { "ID": 4419016,
 "NAME": "", "CAD_NUM": "05:22:000024" }, "geometry": { "type":
 "Polygon", "coordinates": [ [ [...

In PostgreSQL column (cad_num) table value is:

  05:22:000024

Why in the first case does the formatting of the initial value occur? And in the second not?

gdal: release-1900-x64-gdal-3-0-0-mapserver-7-4-0
postgis: postgis-bundle-pg11x64-setup-3.0.0-2


texts: Defining new datum in ArcGIS with seven parameter translation?

texts: Is there a way to define a new datum in ArcGIS with seven parameter translation?


texts: Displaying LiDAR data with too many points in ArcMap?

texts:
I have imported LiDAR data into ArcMap in the shape of a Las DataSet with twelve .las files (contains over 100 million data points).

I believe the file itself is working properly. However, no points are drawn in ArcMap (“Data percentage: 0”) unless I zoom to a very small area. I believe this is a result of the las files containing too much data with over five million points each.

I have tried:

Setting the point limit to five million.
Changing the point density to course.
Filtering out certain points.

But to no avail.

Is there some other way to make the LiDAR data display (at least to some extent) in full view?

Maybe by lowering the number of data points.


texts:
I am new to leaflet.

I am trying to change Leaflet project to EPSG:32624 but it is not working. The map is not showing neither the marker.

var crs = new L.Proj.CRS('EPSG:32634',
    '+proj=utm +zone=34  +datum=WGS84 +units=m +no_defs',
    {
        resolutions: [
            8192, 4096, 2048, 1024, 512, 256, 128,
            64, 32, 16, 8, 4, 2, 1, 0.5
        ],
        //transformation: new L.Transformation(1, XOrigin, -1, YOrigin),
        origin: [0, 0]
        //bounds: L.bounds([218128.7031, 6126002.9379], [1083427.2970, 7692850.9468])
    });
    map = new L.Map('map', {
        crs: crs
    });

L.tileLayer('https://api.geosition.com/tile/osm-bright-3006/{z}/{x}/{y}.png', {
    maxZoom: 14,
    minZoom: 0,
    attribution: 'Map data &copy; OpenStreetMap contributors, Imagery &copy; Kartena'
}).addTo(map);

map.setView([375776.24622049765, 4509140.3932434358], 7);
var m = new L.marker([375776.24622049765, 4509140.3932434358]).addTo(map);
console.log(m);


texts: Importing PostGIS polygon to Esri has negative area and results in null geometry?

texts:
I have data in PostgreSQL with PostGIS that seems fine. I am trying to 'copy' that in to ESRI ArcGIS. I have 4 test rows, 3 import fine. I've also modified a jsfiddle that I found and now has the problem polygon at http://jsfiddle.net/bryanm/nev4hgvo/

fieldnames= ['SHAPE@']
sql = """select st_asbinary(geom),st_astext(geom)
        from asset_utm.dareport_misc
        """
cur.execute(sql)
val = cur.fetchone()
spatRef = arcpy.SpatialReference(4326)
while val != None:
    geom = arcpy.FromWKB(bytearray(val[0]),spatRef)
    geom = arcpy.FromWKT(val[1],spatRef)
    logger.info("Type {0} points {1} parts {2} length{3} area {4}"
                .format(geom.type, geom.pointCount, geom.partCount, geom.length, geom.area ))
    with arcpy.da.InsertCursor( 'VEGDBO.TRAIL_PY',fieldnames) as iCur:
        iCur.insertRow([geom])
    val = cur.fetchone()

The results from the logging are:

2017-07-05 12:24:40,173 - INFO - Type polygon points 20 parts 1 length0.00480354599632 area 1.30829909172e-06
2017-07-05 12:24:40,173 - INFO - Type polygon points 12 parts 1 length0.00575032905258 area 2.10941832519e-06
2017-07-05 12:24:40,188 - INFO - Type polygon points 114 parts 1 length0.00965835903742 area 5.79334045003e-06
2017-07-05 12:24:40,188 - INFO - Type polygon points 58 parts 0 length0.0137069668581 area -9.77088097695e-06

The last row with the -9x10^-6 area is null in the oracle sde based feature. I copied the table into a file geodatabase then ran check geometry, which said it was null.

Other tips I tried was buffering the values by 1 and 0 but that seemed to make things worse - the other geometries became null. I also tried st_simplify with a really small number (0.0000001) but that didn't resolve the issue either - with a larger number I just ended up with 4 points.

First of all, I don't currently know that it is converted to null when I insert it. I've googled for a way to find in arcpy that it is not going to happy.

Is there a better way to check for this than checking for negative area?

Secondly, any idea why this is happening, and how can I resolve it?


texts: Setting up where_clause to select counties one at a time from feature class using ArcPy?

texts:
I need to select a single county at the time from the "counties" feature class. I am using a make feature layer with a where clause.

How do I set up the where_clause to select each of the counties in the row? The field list to be used by the update cursor is "COUNTY". In this list are all the 64 counties names. Here is what I have so far:

# where clause
COname = ??????
where_clause = "COUNTY = '" + COname + "'"
print where_clause
arcpy.MakeFeatureLayer_management (CountyPath, "county_lyr", where_clause)


texts: Can I set the exact height and weight of a letter/character?

texts:
I would like to set the text to have letters 1,7 m high and 1,5 m wide.
Can I do that in QGIS?


texts: Programmatically get selected layers from the QGIS legend?

texts: My QGIS plugin exports all layers from QGIS to KML format. Now I'd like to export only selected layers. However, I can't find any hints in QGIS API. Is it possible to iterate over selected layers in the QGIS legend?


texts:
I've taken multiple gps waypoints and I'm hoping to put them into some sort of GIS program (atm I have GQIS as that was recommended to me) to look at them in map form, but how do I get them there?

I am using GDA94 points and a Garmin GPSMAP 64s.

I have no clue what I'm doing.


texts: Different result with Query builder and the processing 'extract by expression'

texts:
I have a attribute table (xlsx sheet imported with the plugin 'spreadsheet layer'). For example, it contain only 15 lines

I have filtred the data with this request that return 6 lines (correct result) :

  "VALID" != '70' AND "MILIEU"  IN ( 'Entrées de grotte' , 'Grottes et
  cavernes obscures'  , 'Cavité artificielle avec influences humaines' )
  AND  "STRUCTURE" NOT IN ('tunnel','égoûts, canalisation')

Now, if I try to extact data with the same request in the processing modeler, I don't obtain the same result. All data with a NULL value in the field 'Structure' have been filtered. Any idea?

I don't understand why the first code works with the query and not with the extract processing?


texts:
I am trying to complete a select by location, and then Invert my selection, the code runs and does not error, but my selection is not inverted. Not sure what could cause this.

Parks=r"Path"
outputParks=r"Path"
arcpy.MakeFeatureLayer_management(Parks, "Parks")
arcpy.MakeFeatureLayer_management(outputParks, "outputParks")

arcpy.SelectLayerbyLocation_management("Parks", "ARE_IDENTICAL_TO", outputParks, selection_type="NEW_SELECTION", invert_spatial_relationship="INVERT")


texts:
Is it possible to change the precision of a field in the attribute table in ArcMap?

I have a shapefile with X and y coordinates of pointlocations. I would like to change the precision. Is this possible?


texts: Performance Gain Through GIST Index for Point in Polygon Query

texts:
I have two tables: locations(id,region_id,the_geom) and regions(id,the_geom). For each location point I want to determine the region it is located in:

UPDATE locations SET region_id =
 (SELECT id FROM regions
  WHERE ST_Within(locations.the_geom,regions.the_geom)
 );

Does it make sense to build a GIST index on the location points? I'll build an index on the region polygons but I'm unsure about the points. Would it speed up the query?


texts: How to Add Oracle Spatial Tables to QGIS?

texts:
I've been playing qith QGIS Chugiak (2.4.0). I've been attempting to try to add Oracle Spatial tables (point, line and polygon feature classes) to desktop. The same tables I use in ArcGIS. I've managed to get the connection working using the instructions here but 2 things are happening:

Not all tables in my database get populated
When I expand out the the owner groups to get the list of tables, the tables are all greyed out.

I'm not sure why either of these things are happening. In the first issue, I can't see it being a permissions issue as I am logging into the Oracle connection with my Oracle credentials which is already tied to my Oracle account and thus, should grant table access so I'm not sure what's going on there.

In the second issue, I can't figure out why the tables are greyed out. I've tried setting a spatial reference under the SRID field, but that isn't working. Any suggestions?

EDIT:

I got this issue of the greyed out feature tables fixed. I set the SRID field to a relevant spatial reference and also, the 'Primary key column' needs to be populated as well. I used OBJECTID as that is the primary field in my feature classes. So, now I can add my feature classes. I am still having issues getting all of the tables populated. I have several hundred and only 7 of them are loading.


texts: QGIS create raster from binary point layer

texts: I have a point layer with one of the attributes being binary (yes/no).  I want to create a coarse raster which uses the PROPORTION of points within the raster cell which are yes.  So not a count but a count (yes) divided by a count (yes + no). The resolution has to be fairly course because my points are quite spaced - see image.  How can I do this?


texts: Evaluating the geocode precision of US addresses

texts:
I am evaluating various services that take an address in the United States as an input and returns latitude and longitude of the address. Is there a way I can evaluate which service returns the most accurate coordinates?

For example:

Input: 350 Fifth Avenue New York, NY 10118

Output:

Service a: 40.74838, -73.98468

Service b: 40.7483, -73.9853

Is there a way to evaluate which service returned a LAT/LONG closest to the rooftop of the address without manually looking at each set of coordinates on a map?


texts: Errors when sampling Large SpatialLinesDataFrame and sf road layer in R

texts:
I am trying to generate an object that contains regularly sampled points along each road in my road layer (so that there is 0.1 mile/16 meters between each point).  While I can get code to do this on a test road layer I can't get it to work on my curated/edited road layer.

I am trying to figure out what is happening with my road layer and how I can modify it so that regular point sampling operates.  Since this is somewhat specific to my road layer I am not sure if it is appropriate to upload my actual road layer - I am happy to do this if desirable.  Perhaps it is possible to point me in the right direction given the errors?

I have tried, for example, extracting all the line features and saving them in a new object but I don't think I am doing this correctly.  I will not include that code for simplicity and because it may not be germane to my problem here.

Here is the code I am using and the resulting errors:

#Information on road layer:
road

class       : SpatialLinesDataFrame
features    : 854
extent      : 309116.2, 363519, 3767138, 3782434  (xmin, xmax, ymin, ymax)
crs         : +proj=utm +zone=11 +ellps=GRS80   +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
variables   : 2
names       :             rd.name,       rd.type
min values  : 23 NA secondaryRdA3,         4wdA5
max values  : Zuniga Rd localRdA4, secondaryRdA3

class(road)
1 "SpatialLinesDataFrame"
attr(,"package")
1 "sp"

#convert road into a sf object
road.sf  <- st_as_sf(road)
road.sf
Simple feature collection with 854 features and 2 fields
geometry type:  MULTILINESTRING
dimension:      XY
bbox:           xmin: 309116.2 ymin: 3767138 xmax: 363519 ymax: 3782434
CRS:            +proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
First 10 features:
                     rd.name       rd.type                       geometry
1        23 NA secondaryRdA3 secondaryRdA3 MULTILINESTRING ((328799.2 …
2               4wd Rd 4wdA5         4wdA5 MULTILINESTRING ((316419.3 …
3       Abbadie Ln localRdA4     localRdA4 MULTILINESTRING ((349982.3 …
4    Abbotsbury St localRdA4     localRdA4 MULTILINESTRING ((327257.3 …
5  Adamson Flat NA localRdA4     localRdA4 MULTILINESTRING ((342581.2 …
6   Adamsville Ave localRdA4     localRdA4 MULTILINESTRING ((351236.6 …
7     Agua Azul NA localRdA4     localRdA4 MULTILINESTRING ((355929.2 …
8          Alan Dr localRdA4     localRdA4 MULTILINESTRING ((336415.2 …
9         Alder Pl localRdA4     localRdA4 MULTILINESTRING ((347525.1 …
10       Aline Way localRdA4     localRdA4 MULTILINESTRING ((344484.3 …

class(road.sf)
1 "sf"         "data.frame"

################### FAILED CODE ##################################
#here are the errors I am getting from the various functions I have tried:

spsample(road, 16, type="regular")
NULL

spatialEco::sample.line(road, d=16, type="regular", longlat=F)
Error in sp::SpatialPointsDataFrame(lsamp, data = data.frame(LID = rep(as.numeric(lids1),  :
  row.names of data and dimnames of coords do not match
In addition: Warning messages:
1: In mtch != 1:nrow(data) :
  longer object length is not a multiple of shorter object length
2: In sp::SpatialPointsDataFrame(lsamp, data = data.frame(LID = rep(as.numeric(lids1),  :
  forming a SpatialPointsDataFrame based on maching IDs, not on record order. Use match.ID = FALSE to   match on record order

st_line_sample(road.sf, 16)
Error in st_line_sample(road.sf, 16) :
  inherits(x, "sfc_LINESTRING") is not TRUE

st_sample(road.sf, size=16, type="regular")
Geometry set for 5886 features  (with 5870 geometries empty)
geometry type:  MULTIPOINT
dimension:      XY
bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
CRS:            +proj=utm +zone=11 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
First 5 geometries:
MULTIPOINT EMPTY
MULTIPOINT EMPTY
MULTIPOINT EMPTY
MULTIPOINT EMPTY
MULTIPOINT EMPTY


texts:
I'm getting an error while importing metadata to geonetwork (picture below):

TransformerConfigurationException Failed to compile stylesheet. 1 error detected.

I would like to see the details about the error, but I can't find the error log. GeoNetwork is installed on localhost (Tomcat).


texts: Batch load multiple shapefiles to Postgis

texts:
According to the shp2pgsql help, the options drop, append, create and prepare are mutually exclusive. So, if I want to create a table from a shape and then append multiple other shapefiles, I do something like the following, keeping a counter to indicate whether we are in create or append mode.

cnt=0
for shp in $(ls *.shp); do

if [ $cnt -eq  0 ]  ; then

   shp2pgsql -s 27700 -c $shp schema.table_name | psql -h localhost db

else

   shp2pgsql -s 27700 -a $shp schema.table_name | psql -h localhost db

fi
((cnt++))
done

This works, as expected, but I have often wondered if there is a simpler way?


texts: Creating QGIS layers in python console vs stand-alone application

texts:
I am having trouble creating a QGIS vector layer in a stand-alone python script while exactly the same code works fine in python console inside QGIS.  The following code results in a proper displayable layer when executed in the python console:

layer = QgsVectorLayer('/abs/path/buses.shp', 'buses', 'ogr')

However, the same line of code results in an invalid layer in a standalone app (other QGIS functionality seems to work fine):

# PYTHONPATH=/Applications/QGIS.app/Contents/Resources/python
from qgis.core import *

QgsApplication.setPrefixPath("/Applications/QGIS.app/Contents/Resources/python", True)
QgsApplication.initQgis()

layer = QgsVectorLayer('/abs/path/buses.shp', 'buses', 'ogr')
if not layer.isValid():
    print "Layer failed to load: ", layer.lastError() # returns empty string

It seems to me that there is some problem with QGIS data providers in the stand-alone app.  Is there a way to find out what exactly happens inside QgsVectorLayer() or get a list of active QGIS data providers?

I am using QGIS 1.8 on OSX 10.8.3 with python2.7 from macports.


texts: How to remove cache from GWC for a specific layer

texts:
I am trying to remove the cache from a specific layer using curl
so far i got to

curl -v -u user:password -XPOST -H "Content-type: text/xml" -d "<truncateLayer><layerName>workspace:layername</layerName></truncateLayer>"  "http://hosttogeoserver/geoserver/gwc/rest/masstruncate"

I get :  < HTTP/1.1 405 Method Not Allowed
< Date: Fri, 13 Nov 2015 15:37:16 GMT
< Server: Noelios-Restlet-Engine/1.0..8
< Allow:
< Connection: close
< Transfer-Encoding: chunked
< Content-Type: text/plain; charset=UTF-8

Whereas if I do a simple GET request : curl -v -u user:password-H "Accept: text/xml" http://host/geoserver/rest/workspaces/WORKSPACENAME/datastores/DB_NAME/featuretypes?list=available

The response status 200 and outputs the correct information.

I've removed proxy/firewall, etc...

Any ideas?


texts: Converting xyz batch to ESRI grid format using XYZ2GRID 2.1 - batch problems

texts:
I am trying to use Min-Lang Huang's XYZ2GRID2.1-Batch program downloaded from ArcScripts and I keep getting the error

C:\users\public...........xyz Error-Object required

Does anyone know what the problem is? I have created an empty ASCII folder and an empty GRID folder and all the rest of it but not sure what the problem is. Does this script need to sit anywhere special? How does it know that my spatial extension is activated?

I am using ArcGIS Desktop 9.3.1 and Windows Vista.


texts: Choosing the right geo reference system - NorthEast Atlantic

texts:
I am currently trying to calculate fish stock density for various stocks in the NE Atlantic. My habitat areas are too large to use UTM and I was wondering how I would find the best geo-coordinate system to use? My area covers the entire FAO fishing area 27. Any suggestions would be really helpful.

I have changed the meridian as suggested see screenshot. Does this mean that I simply need to change the '10' figure in order to define any new longitude as the central meridian? In other words, currently the central meridian is '10' but if i wanted it to be 15, I simply have to edit this in?

Lambert equal area projection

Saving new CRS projection

Meta Data tab of warped tiff:

Raster --> Projections --> Warp


texts: Features in commercial GIS software missing in FOSS GIS?

texts:
For example, a program I use (TransCAD) has the ability to do travel time matrices based on minimal data input, which is really useful if you're looking at doing back of the envelope type calculations. Another feature it has is the ability to do transit routing and GHG emissions modelling.

What are some of the features you use in non-free/open-source software that you'd like to see implemented in your favorite FOSS GIS software?


texts: QGIS Python SyntaxError file=sys.stderr) File "C:\OSGeo4W64\apps\Python37\lib\site.py", line 177

texts:
I am executing a .py file from a windows .bat file for doing a QGIS Desktop 3.8.3 geoprocessing task.

I am using the following .bat file :

@echo off
call C:\OSGeo4W64\bin\o4w_env.bat
REM Change OSGEO4W_ROOT to point to the base install folder
set QGIS_PREFIX_PATH=C:\OSGeo4W64\apps\qgis
REM Gdal Setup
set GDAL_DATA=C:\OSGeo4W64\share\gdal\
REM Python Setup
set PATH=C:\OSGeo4W64\bin;C:\OSGeo4W64\apps\qgis\bin;%PATH%
SET PYTHONHOME=C:\OSGeo4W64\apps\Python37
set PYTHONPATH=C:\OSGeo4W64\apps\qgis\python;%PYTHONPATH%

REM Launch python job
python C:/Task/Task.py
pause

When I execute this Task.bat in C:/Task/, I am getting the following error.

  File "C:\OSGeo4W64\apps\Python37\lib\site.py", line 177
    file=sys.stderr)
        ^
SyntaxError: invalid syntax
Press any key to continue . . .

There are 2 python files namely Python27 and Python37 in C:\OSGeo4W64\apps\ .

The python-core.bat file in C:\OSGeo4W64\etc\ini\ is :

SET PYTHONHOME=%OSGEO4W_ROOT%\apps\Python37
PATH %OSGEO4W_ROOT%\apps\Python37\Scripts;%PATH%

When I open OSGeo4W.bat in C:\OSGeo4W64\ and type 'python', it is showing me SyntaxError and when I type 'python3', it is showing the version of python as given below.

run o-help for a list of available commands
C:\OSGeo4W64>python
  File "C:\OSGEO4~1\apps\Python37\lib\site.py", line 177
    file=sys.stderr)
        ^
SyntaxError: invalid syntax

C:\OSGeo4W64>python3
Python 3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit
(AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>>

My python code C:/Task/Task.py starts as given below :

from osgeo import ogr, gdal
from gdalconst import *
from qgis.core import *
import qgis.utils, sys
from qgis.gui import *
from PyQt5.QtGui import *
from PyQt5.QtCore import Qt
from PyQt5.QtCore import QFileInfo
app = QApplication([])
QgsApplication.setPrefixPath("C:/OSGeo4W64/apps/qgis", True)
QgsApplication.initQgis()
sys.path.append('C:/OSGeo4W64/apps/qgis/python/plugins')
import processing
from processing.core.Processing import Processing
Processing.initialize()
from processing.tools import *

How can I solve this SyntaxError for executing Task.bat in C:/Task/ ?


texts: How to find the polyline's ID connecting two points in ArcGIS 10.1?

texts:
I have an ArcMap (*.mxd) file with several layers including a point layer and a polyline layer.

The points are connected by polylines in the polyline layer, and I wish to know how to list the connectivity  relation in a table?

For example, if point 1 and 2 are connected by polyline 105, then the relevant row in the resulting table looks something like "1    2    105".


texts:
I am trying to generate a raster with many pixels (10240x10240) using the gdal_grid tool but I have problems choosing the right interpolation algorithm.

I have a dataset with varying point density (ie some places have a lot of points, others have fewer).

With the invdist algorithm, I have "false" results because points far away from a pixel can influence it too much due to their greastest density
I tried to use the radius1, radius2 and max_points parameters, but it seems gdal does not take the closest ones, so the result in meaningless (for my case)
Using the invdistnn algorithm seems to be the solution, but there is no smoothing paramter here which I need because points have a non negligible dispertion.

(I know I could increase the data density artificially where points are missing, but then the computation time becomes very high so this is a solution I would not like to go with.)

Is there a solution to use a smoothing with the invdistnn algorithm ?

I am using GDAL 2.1


texts: Traverse all roads in a city or bounded region

texts:
Is it possible to provide a starting point(lat/long) to a drone and then simulate it to traverse through all the roads in a bounded region.

What I currently have(using OpenStreetMaps) is a list of all intersection points of all the roads in a bounded region(eg: a small sector in a city),

What I want is to make my drone to move through all the roads in that sector automatically(no manual controls) using the lat-long coordinates starting from a particular location, taking minimum possible time.

I have around 5,000 lat-long coordinates(intersection points of roads) in that sector.

Google Direction API provides the use of way-points but only 23 points can be inserted at a time. Besides, Google Maps API require 2 lat-long coordinates(start-point and end-point), what I want is to provide only a starting location and make the drone move along the roads based on some Algo in the backend that covers all roads.

Found a related question here Google Maps City Crawler but the link in solution does not work.


texts: Raster in TOC to geotiff on disk in ArcPy

texts:
I am attempting to automate some ArcGIS10 work w/ netcdf files. Basically, I am able to automate the creation of the netcdf into a raster with arcpy.MakeNetCDFRasterLayer_md. However, I would simply like to then in the same python script save this layer out to a geotiff on my disk so that I can use it later in some analysis. Manually, I know I can go to the TOC right-click and export to geotiff, but it would be great to be able to script this.

I have pasted my arcpy code below for anyone’s reference which works to add the netcdf as a raster to the TOC, but I am stuck after that.

import arcpy

arcpy.env.workspace = "C:/data/test/ArcGISScripting/ArcPy"

mxd = arcpy.mapping.MapDocument("CURRENT")

# Local variables:
crcm_ccsm_cdd_mean_nc = "C:\\Data\\test\\output_files\\crcm_ccsm_cdd_mean.nc"
outputRaster = "C:\\Data\\test\\ArcGISScripting\\ArcPy\\result1_Layer3.tif"

# Process: Make NetCDF Raster Layer
test = arcpy.MakeNetCDFRasterLayer_md(crcm_ccsm_cdd_mean_nc, "result1", "xc", "yc", outputRaster, "", "", "BY_VALUE")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
addLayer = arcpy.mapping.Layer(r"C:\Data\test\ArcGISScripting\ArcPy\result1_Layer3.tif")
arcpy.mapping.AddLayer(df, addLayer, "BOTTOM")
lyr = arcpy.mapping.ListLayers(mxd, "", df)[0].name
#arcpy.mapping.ExportToTIFF(lyr, r"C:\Data\test\ArcGISScripting\ArcPy\result1_Layer3.tif"
del mxd, addLayer
#arcpy.SaveToLayerFile_management("C:\\Data\\test\\ArcGISScripting\\ArcPy\\result1_Layer3.tif", "C:\\Data\\test\\ArcGISScripting\\ArcPy\\result1_Layer3.tif", "ABSOLUTE")
arcpy.RefreshTOC()


texts:
Background:

I want to track the movement of icebergs in ArcGIS. In my FileGeodatabase I have a single Iceberg layer. What I have thought of is to create an ID field for each object; each ArcGIS object that's the same iceberg will have the same ID. They also have a Date field to indicate when they've been observed.

Question:

I would like to only show the latest observation of each object (in ArcGIS terms, for each ID, show the OBJECTID with the latest date).

Using definition query, I've set the where clause to:

SELECT * FROM Icebergs WHERE:
Date = (SELECT MAX(Date) FROM Icebergs m WHERE m.ID = ID)

However no results are returned. Clicking verify shows that the query is syntactically correct. Is this subquery too complex for ArcGIS? Is Definition Query is the correct approach? If not, please point me in the right direction. Later I plan to do this using the Esri JS API, but am using ArcMap to debug the query itself.


texts: Calculating Geoid Height (e.g. EGM2008) for Lat/Lon in PostGIS

texts:
Is there an existing method or best practice for calculating Geoid height for a lat/lon location in PostGIS?  I don't believe there is an "ST_" function.

I know that PostGIS makes use of GeographicLib for some calculations already and that GeographicLib has methods for calculating Geoid height.


texts:
How do I add or join a .lyr to a shapefile?

I have a shapefile of a city's zoning districts and, separately, the .lyr, which shows the correct colors for each zoning classification.

Based on another thread, I tried going to Properties _ Source and changing the Source to point to the shapefile, but then the color distinctions disappeared.

I am new to GIS.


texts: Gdal does not honor proj.4 scaling factor "+k_0"

texts:
My problem:

I'm currently trying to use gdal to write GeoTIFF files. I'm however having problems with using a proj4 string to describe the projection. ("+proj=lcc +k_0=0.001 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +no_defs"). Essentially the scaling factor "+k_0" is ignored by gdal.
I know that instead of +k_0, I can use dataset->SetGeoTransform(adfGeoTransform) (in C++) and set this value manually. However I am surprised that doing this via the proj4 string does not seem to be possible.

If you try the following in python:

Example:

from gdal import osr
p = osr.SpatialReference()
p.ImportFromProj4("+proj=lcc +k_0=0.001 +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +ellps=bessel +x_0=400 +y_0=400 +no_defs")
print p.ExportToProj4()

the output simply drops the +k_0 parameter:

+proj=lcc +lat_1=46 +lat_2=49 +lat_0=47.5 +lon_0=13.333333 +x_0=400 +y_0=400 +ellps=bessel +units=m +no_defs

This results in completely different projected coordinates however and makes it difficult to rely onto.


texts: Prior to image classification, collection of training sites and obtaining a signature file is the most critical step. Field visits with a reliable GPS Unit is usually the first step followed by overlaying the points on the to-be-classified image. Other than this when identifying training sites on satellite images, especially the coarse resolution ones what are the rules normally followed by professionals - like minimum number of pixels in a training site, distance between two training sites for same/different classes, how improvements are made after preliminary classification etc. Standard references are also welcome.


texts:
I need to optimize the existing route in my city for my undergrad thesis. So the idea I had is to use a tool to analyze the current existing route's accessibility and comparing it to the accessibility of the new route I'm about to design by trial and error (Designing multiple new routes and finding the one with the biggest accessibility differences). But I'm completely blind for this tool.
Is there any existing tool I can use? Or a plugin? Or do I need to create my own algorithm in the create new script menu?


texts: Workaround for AxMapControl painting issue (only problematic at certain screen resolutions)?

texts:
I have a .NET application dependent on ArcGIS 9.3.1 which I am upgrading to ArcGIS 10.

I have an AxMapControl sat in a panel in a split container on a tab. The Dockstyle is set to Fill.

Since upgrading to ArcGIS 10 from 9.3.1, the painting of the map inside the map control has become inconsistent.

The initial size of the map as painted is 'wrong' and the map ends up only filling, say, half the panel, and a second vertical scrollbar appears beside the map. This second scrollbar is a glitch; you cannot use it, and it is eventually painted over. The issue with the the map seems to be painting only: all other information presented (for example in the legend) is as if the map were painted correctly across the whole panel.

After a few operations (such as zooming in/out etc.) the map painting seems to 'catch up' and the map is painted correctly, filling the whole panel. This 'catching up' seems to get faster each time I run the application which leads me to suspect some sort of timing or caching issue.

There is no consistency about when the problem manifests. Sometimes adjusting the screen resolution affects it, other times not.

This seems to be a problem internal to the control. Does anyone have a workaround (other than reverting to the 9.3.1 DLLS)?

UPDATE: I have tried the following without success:

Setting the dock property of the mapcontrol to fill (http://forums.esri.com/Thread.asp?c=159&f=1707&t=209034#631591)
Manually resizing in the resize event of the mapcontrol
Deactivating and reactivating the ActiveView
Swapping the map layers around
Forcing lots of Full Refreshes
Invalidating the map


texts: Calculate Field throws an error when parsing python with backslashes

texts:
I have an FTP I need to access through a python script.  It works fine with forward slashes, but I don't need it to be read as an html. Below are both sets of code:

Working:

arcpy.CalculateField_management(outShp, "link", '"//ftp/raw" + time.strftime("%Y%m%d") + "/files/boundary" + !NAME! + ".shp"'

Invalid with backslashes:

arcpy.CalculateField_management(outShp, "link", '"\\\\ftp\\raw" + time.strftime(%Y%m%d") + "\\files\\boundary" + !NAME! + ".shp"'

Pretty new, but seems like it should work with the escapes, and does when used in the Calculate Field tool - just doesn't translate to my script.


texts: Sources for long-term NDVI datasets with coverage to west Asia?

texts: I'm in need for long-term NDVI dataset for my research in bio-climatology, looks like GLCF GIMMS are no longer available online, any suggestions/workarounds to get the data.


texts:
I have a table that has its geometry type as geometry and I need to convert it to Polygon for ogr2ogr to read it. Can't figure this out.

I have tried

ALTER TABLE oak_all_buffers ALTER COLUMN geom type geometry(MultiPolygon, 102689);

error

    ERROR:  Geometry type (MultiPolygon) does not match column type (Polygon)
    ********** Error **********

    ERROR: Geometry type (MultiPolygon) does not match column type (Polygon)
    SQL state: 22023

and

ALTER TABLE oak_all_buffers ALTER COLUMN geom type geometry(Polygon, 102689);

ERROR:  Geometry type (Polygon) does not match column type (MultiPolygon)
********** Error **********

ERROR: Geometry type (Polygon) does not match column type (MultiPolygon)
SQL state: 22023


texts: How to display map tips only for one layer?

texts:
I'm relatively new to QGIS.

I added the hover labels to my map using the Map Tips. My issue is that it makes hover labels for EVERY layer on my map. I only want it to make labels for one layer.

Is there a way to fix this like assigning Map Tips to only apply to one layer?


texts:
I wanna to contribute for algorithms, that's why i chose pgrouting can anyone help me out and suggest me more organisations to which i can contribute.
Thnks.


texts: CartoDB: Import Error: host and port must be set

texts:
I have setup cartodb in one of the server (CentOS). Everything worked fine in development environment, but when I move to production environment I couldn't import any file. Getting the below error in production.log

Ignoring statsd, because there were a error loading the config
[table#create error] === #<RuntimeError: host and port must be set>

Below the importer job log,

2014-06-17T10:44:54Z || Before dispatch
2014-06-17T10:44:54Z || new_importer()
2014-06-17T10:44:54Z || Fetching datasource public_url metadata for item id http://cartodb.s3.amazonaws.com/static/TM_WORLD_BORDERS_SIMPL-0.3.zip
2014-06-17T10:44:55Z || File will be downloaded from http://cartodb.s3.amazonaws.com/static/TM_WORLD_BORDERS_SIMPL-0.3.zip
2014-06-17T10:44:55Z || Before run.
2014-06-17T10:44:57Z || Starting import for /tmp/imports/20140617-14039-17m47zj/TM_WORLD_BORDERS_SIMPL-0.3.zip
2014-06-17T10:44:57Z || Unpacking /tmp/imports/20140617-14039-17m47zj/TM_WORLD_BORDERS_SIMPL-0.3.zip
2014-06-17T10:44:57Z || Importing data from /tmp/imports/20140617-14039-166nffe/tm_world_borders_simpl-0.3.shp
2014-06-17T10:44:58Z || Detected encoding WIN1252
2014-06-17T10:44:58Z || Using database connection with {:adapter=>"postgres", :encoding=>"unicode", :host=>"localhost", :port=>5432, :database=>"cartodb_user_d13717ee-414f-4c29-a4e9-19e1d48c197b_db", :username=>"postgres", :user=>"cartodb_user_d13717ee-414f-4c29-a4e9-19e1d48c197b"}
2014-06-17T10:44:59Z || ogr2ogr call:      OSM_USE_CUSTOM_INDEXING=NO PG_USE_COPY=YES PGCLIENTENCODING=WIN1252 SHAPE_ENCODING=WIN1252 /usr/local/bin/ogr2ogr -f PostgreSQL PG:"host=localhost port=5432 user=cartodb_user_d13717ee-414f-4c29-a4e9-19e1d48c197b dbname=cartodb_user_d13717ee-414f-4c29-a4e9-19e1d48c197b_db password=71fc094fbf3ac1d38de9cc90ac3eeed821df33f0cartodb_user_d13717ee-414f-4c29-a4e9-19e1d48c197b active_schema=cdb_importer" -t_srs EPSG:4326  -lco DIM=2 -lco PRECISION=NO /tmp/imports/20140617-14039-166nffe/tm_world_borders_simpl-0.3.shp  -nln cdb_importer.importer_6cd598c4f60c11e383fd00185162f697 -nlt PROMOTE_TO_MULTI
2014-06-17T10:44:59Z || ogr2ogr output:
2014-06-17T10:44:59Z || ogr2ogr exit code: 0
2014-06-17T10:44:59Z || Georeferencing...
2014-06-17T10:44:59Z || Disabling autovacuum for "cdb_importer"."importer_6cd598c4f60c11e383fd00185162f697"
2014-06-17T10:44:59Z || Creating the_geom from wkb_geometry column
2014-06-17T10:44:59Z || replace empty strings by nulls?
2014-06-17T10:44:59Z || no string column found, nothing replaced
2014-06-17T10:44:59Z || geometrifying
2014-06-17T10:44:59Z || casting wkb_geometry to geometry
2014-06-17T10:44:59Z || Converting to 2D point
2014-06-17T10:44:59Z || Total processed: 246
2014-06-17T10:44:59Z || Total processed: 246
2014-06-17T10:45:00Z || FINISHED: Converting to 2D point
2014-06-17T10:45:00Z || geometrified
2014-06-17T10:45:00Z || Renaming column wkb_geometry TO the_geom
2014-06-17T10:45:00Z || Enabling autovacuum for "cdb_importer"."importer_6cd598c4f60c11e383fd00185162f697"
2014-06-17T10:45:00Z || Georeferenced
2014-06-17T10:45:00Z || Finished importing data from /tmp/imports/20140617-14039-166nffe/tm_world_borders_simpl-0.3.shp
2014-06-17T10:45:00Z || Proceeding to register
2014-06-17T10:45:00Z || Before renaming
2014-06-17T10:45:00Z || Before moving schema
2014-06-17T10:45:00Z || Before persisting metadata
2014-06-17T10:45:01Z || Dropping table tm_world_borders_simpl_0_3
2014-06-17T10:45:01Z || Import Error: host and port must be set
2014-06-17T10:45:01Z || After run.
2014-06-17T10:45:01Z || After dispatch
2014-06-17T10:45:01Z || SUCCESS!

Which host/port are missing and where to set?


texts: Adding a google Maps Layer to OpenLayers 2.14

texts:
I have a well-working map consisting of several layers in OpenLayers retrieved from geoserver. Now, I want to add a Layer from google maps.
Since I did not find any working explaination in the web - how do I achieve this?


texts: Cartesian / Meters issue while making ellipses

texts:
I am trying to create an ellipse to be stored in my django project. I'm provided a lat, lon, major/minor in Nautical Miles, and an orientation in degrees. Between this post and the docs I've been able to piece together the below code, however, when I compare the final outputs from an export of my django to what I have generated from the arcgis "Table to Ellipse tool" they don't match. I've tried to account for Arc's different methodology of using azimuth and clockwise rotation in my code but the orientation and scale don't match. I think the issue stems from starting in cartesian space as meters, converting to WGS-84, and then translating origin position. Here is my code with random test data.

from functools import partial
from shapely.ops import transform
from shapely.affinity import translate
from django.contrib.gis.geos import GEOSGeometry
from django.contrib.gis.geos import MultiPolygon
import pyproj

lat = -78.150511
lon = 39.190978
major = 13.114 #Nmi
minor = 8.58 #Nmi
metersconversion = 1852
smajormeters = (major * metersconversion) / 2
sminormeters = (minor * metersconversion) / 2
orient = 61.9 #Degrees Rotation
#Ellipse parameters
ellipse = ((0, 0), (smajormeters, sminormeters), orient)
#create circle with radius of 1
circ = shapely.geometry.Point(ellipse[0]).buffer(1)
#create ellipse along x and y
ell = shapely.affinity.scale(circ, int(ellipse[1][0]), int(ellipse[1][1]))
#Rotate ellipse
elrv = shapely.affinity.rotate(ell, 90 + (ellipse[2] * -1))
project = partial(pyproj.transform, pyproj.Proj(init='epsg:3857'),
pyproj.Proj(init='epsg:4326'))
geom4326 = transform(project, elrv)
transgeom = shapely.affinity.translate(geom4326, lon, lat)
ellpoly = GEOSGeometry(transgeom.wkt)
#Define the srid/epsg as 4326 for ingest into django postgresql
ellpoly.srid = 4326

What has gone wrong?


texts: SAS Planet cached tiles for use in Leaflet map

texts:
I downloaded some tiles from ArcGIS.Imagery with SAS Planet
The tiles is in folders like

ArcGIS.Imagery/z2/0/x1/0/y2.jpg
ArcGIS.Imagery/z13/2/x2624/1/y1750.jpg

And I use this TileLayer to display it on Leaflet.

  L.TileLayer.customTiles = L.TileLayer.extend({
            getTileUrl: function (t) {
                console.log('coords', t);
                z = t.z,
                    x = t.x,
                    y = t.y,
                    z0 = t.z,
                    x0 = Math.floor(t.x / 1024),
                    x1 = Math.floor(t.x % 1024),
                    y0 = Math.floor(t.y / 1024),
                    y1 = Math.floor(t.y % 1024)

                const d = `ArcGIS.Imagery/z${z0}/${x0}/x${x1}/${y0}/y${y1}.jpg`;
                console.log('d', d);
                return d;
            }
        });
        L.tileLayer.customTiles = function () {
            return new L.TileLayer.customTiles();
        }
        L.tileLayer.customTiles().addTo(map);

But it looks like there is wrong coordinates. Here is what I got as a result.

Any ideas how to fix it?


texts:
I can't remember if AutoCAD supports reprojecting layers from a projection/datum to an other projection/datum.

My first question is : could it be done on AutoCAD?

Second question, which matters more to me -since I'm a GIS person - is : could it be done by ArcGIS or QGIS or any of the other GIS softwares ?

I've done some research on the internet, it turns out that a georeferencing of a CAD file is possible in ArcGIS, but what's a georeferencing operation worth if the spatial reference of the file is wrong in the first place ?

Note : of course one can think of using the interoperability tool in ArcGIS, convert a DWG to a SHP, reproject the SHP then finally re-convert to DWG. I'm looking for a simplier solution.

Detailed description : we have a DWG file in the "Lambert Conical Conformal Zone 1 / Clark 1880" for Morocco, and we want to know if it can be possible to reproject it into "UTM Zone 30N" (WGS84). This question would've not had posed any problem if it was a Shapefile, as you would imagine. Thus, the question is in short : what are known ways to reproject DWG files from their spatial reference to a target spatial reference.


texts: How to do a raster reprojection in QGIS with custom map parameter

texts:
I'd need to reproject a raster dataset in a new customized projection, but I'm fairly new at both QGIS and GDAL and I'm not exactly sure about how to proceed and what the exact syntax is.

The starting raster dataset is in EPSG:3035 and I'd like to reproject it in a Lambert Conic Conformal projection of which I know center latitude, center longitude, true latitude 1, true latitude 2, standard longitude and pole latitude and longitude (the latter are standard: 90 and 0 degrees).

I guess that I should use "Raster"-->"Projection"-->"Reproject" under QGIS, but I'm not quite sure on how to express it in proj4 syntax.

Thanks for any help!


texts: How to solve Error 000840 from arcpy.JoinField_management?

texts:
As part of a larger Python script, I am attempting for my last step to join an Excel file to a spatial layer.  I was able to perform this step just fine in ArcGIS Pro Model Builder, but when I attempt to do it with a stand-alone script in Python I get an error (below).
I imported pandas and arcpy at the beginning.  I set the workspace to 'G:\Project.'  Then there are intervening steps with the Excel file that all ran just fine.  But here is the problematic script.
# Join new Excel file to spatial layer
inFeatures = r'G:\Project\spatiallayer.lyrx'
joinTable = r'G:\Project\Excel.xlsx'
joinField = "OBJECTID"

joined_table = arcpy.JoinField_management(inFeatures, joinField, joinTable, joinField)

And here is the error I get:
ExecuteError: Failed to execute. Parameters are not valid.
ERROR 000840: The value is not a Table View.
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
Failed to execute (JoinField).

I originally had a layer package, but Python would not recognize that, so I created a new layer from that package.  And the error above is what I get when I attempt the Join Field with that new layer.  What makes this so baffling is that it works in ModelBuilder using the same inputs.  I don't understand why it will not accept my layer to join with my Excel file.  How can I solve this error?
I am running Python 3.6.1 and ArcGIS Pro 2.6.1.


texts: Filter layer by using symbology as a toggle

texts:
Is there a way to customize ArcGIS Desktop 10.7.1 so that users can filter layers by simply clicking the layer's symbology in the table of contents?

For example, is there a way to add a toggle/checkbox to each level in a layer's symbology that users could click to include/exclude the features from the map?

Or alternatively, select the symbol's label/text to include the features in the map. And deselect the symbol's label/text to exclude the features.
Or another way of approaching it might be: enable an option in the layer's properties that allows users to generate a definition query by simply clicking in the symbology.

This functionality would add efficiency to day-to-day work by eliminating the need to write tedious definition queries (WHERE clauses with IN operator) just to include/exclude a certain type of feature from a layer.


texts: OpenLayers - Add points at vertices of line while drawing line

texts:
I would like to create points while I draw a line on my map in OpenLayers. That is, at each vertices of a line and at the end points there is a shape like a square, circle, or diamond.

I have this partially working. I can add points at the vertices of a line while I'm drawing it but they are one off. The point/shape doesn't render at the beginning of the line but at the second point in the line as you can see below. Also, the end of the line has two points/shapes; you can't see it since they are stacked.

Here is my code that adds the point/shape as I'm drawing my line (see the point callback in my DrawFeature control setup).

var path_control = new OpenLayers.Control.DrawFeature(
     vector_layer,
     OpenLayers.Handler.Path,
     {
          callbacks: {
               point: function( point ) {

                    console.log(point);

                    point_collection.push(point);
                    var feature = new OpenLayers.Feature.Vector( point );

                    vector_layer.addFeatures([feature]);
                    vector_layer.redraw();

               }

          }
     }
);

console.log(point) gets called at all the correct times. When I output my point_collection variable that I setup to figure out what is going on, the points are different than those output by console.log(point) from with the point callback.

I have a feeling that it has something to do with when I'm actually adding the points to the vector layer but can't figure it out.


texts: Opening personal geodatabases (*.mdb) in QGIS 3.4.1?

texts:
I was successfully able to open *.mdbs in my previous version of QGIS on Windows 10 64 bit, but since the Madeira update (3.4.1, 64 bit) they no longer work.

Getting a "Layer is not valid" error with no given reason.

I have added them as custom variables through QGIS:

OGR_SKIP
value: ODBC

PGEO_DRIVER_TEMPLATE
value: DRIVER=Microsoft Access Driver (*.mdb, *.accdb);DBQ=%s

I've tried adding these variables in through Windows advanced systems settings/environments as well, same variables as above.

Gone through multiple restarts and combinations of variable settings.

Microsoft Access database engine (64) is installed.

Is it the update? Am I doing something wrong?

Alternatively, is there another way to access ESRI personal geodatabases through QGIS? I haven't been successful using any other methods I've found so far but any tips would be great!


texts: Simplifying a polygon, result must contain original

texts:
I have some high-res polygons, without enclaves or exclaves that I want to simplify or smooth without loosing any part of the original.
Using NettopologySuite I could do polygon.Buffer(x) and then DouglasPeuckerSimplifier with DistanceTolerance=x
Perf doesn't matter since this is done during build.


texts: Reconciling Features with different field types

texts:
I have two feature classes who mostly have the same fields. However, some of the fields do not share the same field types. For example, one has a created date field that is a Date while the other feature class has a created date field that is a String. This is repeated on several fields with doubles, shorts, etc.. Normally, no problem, I can use a join and the field calculator.

However, these two feature classes represent before and after GPS data collection. So one feature might be a few feet off of it's counterpart. One last constraint is that the info and data must be migrated to one of the feature classes in it's existing database with existing relates. At least I think that is a constraint because the database is setup in a certain way to work with the surveying software.

Data needs migrated to this feature using the field types existing here:

Data from here needs migrated and in some cases converted between field types. Actual location from this one also needs to be used:

Are there any tools or recommendations to reconcile these two feature classes into a single feature class as described? I'm guessing it's going to be a painful multi-step process.


texts: Cropping openstreemap tiles to desired bounding box

texts:
Given the solution for downloading a cluster of tiles in

https://stackoverflow.com/questions/28476117/easy-openstreetmap-tile-displaying-for-python

and

http://wiki.openstreetmap.org/wiki/Slippy_map_tilenames#Python

what is a good way of EXACTLY cropping the resulting image to the DESIRED bounding box using Python?

Is interpolation necessary?

Note that when downloading OSM tiles you usually get returned a bounding box that is larger or equal the desired one due to the discreet nature of the tiles.


texts: How to register as versioned using SDELayer?

texts:
This is the normal process using ArcCatalog :

I want to automate this process using sdeLayer or any tools - does anybody have an idea as to how to do it?

there is an application sdelayer (sdelayer.exe) in "C:\Program Files\ArcGIS\ArcSDE\pgexe\bin" i want to create a command line to register featurelayer

sdelayer -o register -l <table,column> -e <entity_mask> -t <storage_type>
{[-C <row_id_column>[,{SDE|USER}[,<min_ID>]]]} [Spatial_Index]
[-E <{empty | xmin,ymin,xmax,ymax}>] [{-R <SRID> | [Spatial_Ref_Opts]}]
[-S <layer_description_str>] [-k <config_keyword>]
[-i {<service> | <port#> | <direct connection>}] [-s <server_name>]
[-D <database_name>] -u <DB_user_name> [-p <DB_user_password>] [-q]


texts: Failed to update feature through WFS with QGIS

texts:
I am trying to update a polygon shapefile via a WFS service created using GeoServer. The client used for updating WFS is QGIS. So far I have been able to create a new feature successfully but I'm unable to update/delete an existing feature. Every time I try to update the existing feature, a new feature is created and the old feature is also retained. While deleting a feature does not work.

I am unable to figure the reason why.

Note: Editing of WFS has to be through QGIS


texts: How to add a legend, banner to shp file exported in PDF

texts:
I need add a legend, banner and signature to my shp file which is exported in PDF.
My code:
HashMap<String, Object> params = new HashMap<>();
        params.put(ShapefileDataStoreFactory.URLP.key, URLs.fileToUrl(new File("D:/var/shp/poli_cercanos.shp")));
    DataStore ds = DataStoreFinder.getDataStore(params);
    SimpleFeatureCollection fc = ds.getFeatureSource(ds.getTypeNames()[0]).getFeatures();

MapContent mapContent = new MapContent();
mapContent.setTitle("Quickstart");
Style style = SLD.createSimpleStyle(fc.getSchema());
Layer layer = new FeatureLayer(fc, style);
mapContent.addLayer(layer);

ByteArrayOutputStream baos = new ByteArrayOutputStream();

int w = 1000;
ReferencedEnvelope bounds = fc.getBounds();
int h = (int) (w * (bounds.getHeight() / bounds.getWidth()));
BufferedImage bufferedImage = new BufferedImage(w, h, BufferedImage.TYPE_INT_ARGB);
Graphics2D g2d = bufferedImage.createGraphics();

mapContent.getViewport().setMatchingAspectRatio(true);

mapContent.getViewport().setScreenArea(new Rectangle(Math.round(w), Math.round(h)));
mapContent.getViewport().setBounds(bounds);

g2d.setRenderingHint(RenderingHints.KEY_ANTIALIASING, RenderingHints.VALUE_ANTIALIAS_ON);

Rectangle outputArea = new Rectangle(w, h);

GTRenderer renderer = new StreamingRenderer();
LabelCacheImpl labelCache = new LabelCacheImpl();
Map<Object, Object> hints = renderer.getRendererHints();
if (hints == null) {
    hints = new HashMap<>();
}
hints.put(StreamingRenderer.LABEL_CACHE_KEY, labelCache);
renderer.setRendererHints(hints);
renderer.setMapContent(mapContent);
renderer.paint(g2d, outputArea, bounds);
ImageIO.write(bufferedImage, "png", baos);

Document document = new Document(PageSize.A2, 20, 20, 20, 20);
PdfWriter.getInstance(document, new FileOutputStream("D:/var/shp/poli_cercanos.pdf"));
document.open();
Image image = Image.getInstance(baos.toByteArray());
document.add(image);
document.close();`

My pdf generates the following image from the shp file:

Then the following would be to add a legend and a signature to the image, like this:


texts: OpenLayers Popup Anchored Offset Problem

texts:
I have an animation where I move a point along a line. At certain spots of that line I open a popup to show a picture. I want the popup to be centered on the point. It seems to work for some points and not for others. Sometimes the popup is centered and sometimes it is just offset off the point. I tried to set that the relative position is "tr", but that did not help(maybe I did it wrong?):

    var popup_point = new OpenLayers.LonLat(fence_point.x,fence_point.y);
    var anchor = {'size': new OpenLayers.Size(0,0), 'offset': new OpenLayers.Pixel(-100, -100)};

    popup = new OpenLayers.Popup.Anchored("fence",
                       popup_point,
                       new OpenLayers.Size(200,200),
                       contentString,
                       anchor,
                       false);

    // popup.autoSize = true;
    map.addPopup(popup);


texts:
When trying to export a map as an MPK, ArcMap gets stuck at "Preparing Data (Folder Size 0KB)". I have tried a couple different MXDs and they all have this problem.

I am working off a server with a not-so-great connection, but the connection isn't that bad.  Any thoughts on what could be causing this?


texts: Is there a plugin to import places of interest (restaurants, hotels, parks, etc.) into QGIS from Google Maps?

texts: I am developing a workflow to identify important/interesting locations contained within a pre-defined area. I would like to know if anyone can indicate a plugin or method to import locations from Google Maps that I can identify if they are contained within my AOI.


texts: Single sided buffer produces no result in QGIS?

texts:
I've mostly used ArcGIS Desktop but am trying to use QGIS to do a one sided buffer since I only have a basic license.

I have a coastline and want to do an offshore buffer polygon. The line file is projected to a state plane system consistent with my desired buffer units.

I've opened "Single sided buffers (and offset lines) for lines" put in my coastline shapefile line file, selected single side buffer, selected geometry, input buffer distance, selected right, and specified an output shapefile. When I Run this nothing happens...

What am I missing?


texts: Extracting Overlapping Parts of Single Polygon Layer using QGIS?

texts:
I have a single polygon layer with over 200 buffers with a lot of overlapping areas. I want to use QGIS to create a new layer representing the overlap. I've seen several questions that are similar, but haven't been able to figure it out yet. I'm new to QGIS, so I may be overlooking a simple solution.
I'd like to save out the red overlap from this example:

From this question:
Remove overlapping buffers
Similar to this question, but I'm using QGIS 2.8.1
how to extract overlap from non-dissolved buffer zones in qgis 1.7.4


texts:
I have a problem connected with my Raster Dataset projection. Two of its bands cannot be shown in some scale - only under 1:1500, when I use 'Zoom to raster resolution' option.
Histograms for all 3 bands are correct, I even calculated statistics and did pyramids.

Below you can find screens of it (3rd band is presented same as band 2):

I would be grateful if you could help me :)

PS. I store them (R.D.) in a file geodatabase (.gdb).


texts: Able to plot on map but having trouble with st_contains

texts:
I am new to working with GIS, and am building a rails app that uses the rgeo gem and PostGIS.

I am trying to use geospatial data to determine the school district for a house. The house is represented by a point, and the school districts are represented by multipolygons. I have successfully plotted both on a map, and all appears to be working correctly, with the point falling inside of one of the polygons.

However, when I try to run the following code to determine which district the house is in, I get back an empty array:

>>  point = @property.coords
>> ES_Zone.where{ st_contains(geometry, point) }
=> #<ActiveRecord::Relation []>

I think that the SRIDs are the same on the points and polygons, because the factories for each show SRID 3785.

>> point.class
=> RGeo::Geos::CAPIPointImpl
>> point.factory
=> #<RGeo::Geos::CAPIFactory:0x860c5838 srid=3785 bufres=1 flags=8>

>> ES_Zone.first.geometry.class
=> RGeo::Geos::CAPIMultiPolygonImpl
>> ES_Zone.first.geometry.factory
=> #<RGeo::Geos::CAPIFactory:0x84889648 srid=3785 bufres=1 flags=8>

Any ideas on why this is happening, and suggestions for how to fix it?