texts: pgRouting 2.0 slow execution of A star algorithm

texts:
Recently, I wanted to test pgRouting ver. 2.0. I created the database structure with osm2po and imported UK map from geofabrik. The data was imported correctly.

When I execute a sample query in order to find the shortest path with dijkstra algorithm I receive the result in ~70 sec. When I do the same with astar algorigthm I get the result in over 145 sec. Theoritically, the result of astar should be faster than dijkstra's one but it's not.

Queries I execute:

Dijkstra

SELECT seq, id1 AS node, id2 AS edge, cost
    FROM pgr_dijkstra(
        'SELECT id::integer, source::integer, target::integer, cost::double precision AS cost
             FROM hh_2po_4pgr'::text, 60, 80, false, false);

Astar

SELECT seq, id1 AS node, id2 AS edge, cost
FROM pgr_astar(
    'SELECT id AS id, source::integer, target::integer, cost::double precision AS cost, x1, y1, x2, y2
        FROM hh_2po_4pgr', 60, 80, false, false);

I checked the Astar implementation in pgRouting and it seems that it uses this heuristic: (|dx|+|dy|)/2.

Am I doing something wrong in my queries or there is something wrong with pgr_astar() implementation in pgRouting?


texts: How do I format numbers with thousand commas in custom infobox in Carto?

texts:
For example, I'm calling a number into the infobox using the following code:

{{med_value}}

The value is 1155000, but I want it to display as 1,155,000.

Is there a way to do that without having to reformat the column as a string and adding the commas that way?


texts: Merge vector points with point name into field name ArcGIS

texts:
I am merging several points gathered on Google Earth after converting KML into a shapefile. I would like to have in the output file (merged) a field indicating the file name of the corresponding source points (the original single points that I want to merge).

Is that possible in ArcGIS?


texts: arcpy.SelectByAttribute 'Parameters are not valid' error

texts:
I am trying to run a python script to select 5 random points in an attribute table and to loop through every 100 points to do this. When I run the function SelectRandomByCOunt within a python shell, it works perfectly fine but I am trying to embed it in this For loop and it doesnt work. I keep getting this error mesg:

ExecuteError: Failed to execute. Parameters are not valid.
The value cannot be a feature class
ERROR 000840: The value is not a Raster Layer.
ERROR 000840: The value is not a Mosaic Layer.
Failed to execute (SelectLayerByAttribute).

The code

import arcpy
import os
from arcpy import env

# path name of the feature class
ws = r"S:\Venu\APWRAGOEAProject\APWRAGOEA.gdb"
fcname = r"PracticePythonCode"
fc = os.path.join(ws, fcname)

# replace 'xxx' with file path

# get count of features
num = arcpy.GetCount_management(fc)
print("There are "+str(num)+ " rows in your feature class.")

# start counters
count = 0
next100 = 1

with arcpy.da.SearchCursor(fc, ['OBJECTID']) as scursor:  ## column 1 name [0], column 2 name [1]

      for srow in scursor: # iterate through rows

        if count < (next100 * 100):
                      col1 = srow[0] ## grab value from column 1

      # do processing of these 100 records here using col1, col2 variables
                      def SelectRandomByCount (layer, count):
                                import random
                                layerCount = int (arcpy.GetCount_management (layer).getOutput (0))
                                if layerCount < count:
                                    print "input count is greater than layer count"
                                    return
                                oids = [oid for oid, in arcpy.da.SearchCursor (layer, "OID@")]
                                oidFldName = arcpy.Describe (layer).OIDFieldName
                                path = arcpy.Describe (layer).path
                                delimOidFld = arcpy.AddFieldDelimiters (path, oidFldName)
                                randOids = random.sample (oids, count)
                                oidsStr = ", ".join (map (str, randOids))
                                sql = "{0} IN ({1})".format (delimOidFld, oidsStr)
                                arcpy.SelectLayerByAttribute_management (layer, "", sql)

                      SelectRandomByCount (fc,5)

                      # increase count by 1
                      count = count + 1

        else:

              # move to next group of 100
                next100 = next100 + 1

del srow, scursor


texts: Import GRD and count days with temperature below zero, using python

texts: Minimum temperature of a year is given in GRD file. I have to read that file and then calculate the frost day (less than zero degrees celcius) from it in python.


texts: Is it necessary/advisable to uninstall QGIS 1.8 first before installing 2.0?

texts: Who has installed QGIS 2.0 on Mac OSX 10.6.8 (Snow Leopard)? Is it necessary/advisable to uninstall QGIS 1.8 first? Have you seen any problems importing QGIS 1.8 projects to 2.0?


texts: Creating linear connections among points

texts:
I would like to create in an automatic way a linear feature representing the connections among all points using ArcMap or Python.
That is I would like that each point is connected to the others one like the attached pic.

import arcpy
import os
arcpy.env.workspace = r"C:/test/Test.gdb"
fc_list = arcpy.ListFeatureClasses()
array = arcpy.Array()
pnt = arcpy.Point()
for record in fc_list:
    polyline = arcpy.Polyline(pnt)
arcpy.CreateFeatureclass_management("C:/test/Test.gdb", "linea", "POLYLINE")
print("finished")

The script finishes, but the "linea" created is empty. In the "Test.gdb" there is only a point feature class called "gps" which has 17 points that have to be connected each other.


texts: Checking availability of ET Geowizards license

texts:
I would like to run a script that uses some tools from the ET Geowizards toolbox. I have a valid license for it but sometimes the license server at work is off. If the server is off then some of the ET Geowizards tools mentioned above are limited to processing only 100 features.

I can import the ET Geowizards toolbox via arcpy.ImportToolbox with no problem, import works even if the license is not available. Bad thing is that I don´t know if the server is running or not at the time of arcpy.ImportToolbox. I must wait for the script execution if it raise an error or not. If an error has been raised, then the server with the license is off.

Is there some python/arcpy way of checking the availability of the ET Geowizards license on the server? Something like arcpy.CheckExtension for ArcGIS Extensions.

Running on ArcGIS for Desktop 10.1 Basic + SP1, ET Geowizards 10.2 for ArcGIS 10.1, Win8 64 bit


texts: Project Polygon using ArcGIS API for JavaScript

texts:
I recently did a javascript application that has since stopped working, well partially. The users enters an address, then they are presented various voting related information State Rep, State Senator, voting precinct, school dist etc. When the user clicks 'show district', the respective district (State Rep, School Dist etc) polygon should highlight, and indeed it did when deployed.

Just recently I was notified it did not work anymore. It seems to be a problem with the basemap being in one spatial reference and my reference layers being in another. I can possibly get my ref layers changed, but this will require an act of congress. I was trying to reproject my reference polygons via GeometryService to project onto my basemap, but got nowhere. My geometries did not return in the callback.

Current Prod: http://argis.ualr.edu/districtfinder/

Address (copy/paste): 1 Capitol Mall, Little Rock, Arkansas, 72201


texts: Hiding QGIS toolbars through python code?

texts:
for eg.
iface.rasterToolBar().setVisible(False)

this works from python console of QGIS 2.6.1 but does not work from plugin code.


texts:
I have a huge problem that I can't solve even after trying everything I can think of.

I am trying to do some analysis of the Lake Chad basin in QGIS (2.16.3). I downloaded Landsat 8 data from Earth Explorer, however I need 4 different satelite images to cover the whole basin. I did atmospheric correction and created a bandstack, the problem that now appears is that I have a "black border" around the single images.

I tried to fix this (as read in another post in here) by going to layer properties -> transparency -> additional no data value = 0 but it didn't work (I also tried raster -> analysis -> near red), so instead I took the single image, and did raster -> miscellaneous -> merge -> and checked "no data value". This worked to remove the black borders, but now when I try to merge the four single images into one, the black borders appear again.

Does anyone know how to fix this?

I've been trying to do this now for days, but I really don't know how to help myself at this point. I will add some images to clarify what I'm talking about.


texts: OpenLayers fill opacity increases after repositioning

texts:
I have created a test project where i'm trying to learn how to work with OpenLayers. I try to do a get request to a remote geoservice. Every time I reposition the filled polygons on my map get there opacity increased. So I think the polygones are put on top of each other. Here you can view the URL that I from where i get my JSON data. This is done over a Node JS server on localhost beceause the original link doesn't support CORS.

var text =
      "http://localhost:3030/map/https://geoservices.landbouwvlaanderen.be/PUBLIC/wfs?service=WFS&request=GetFeature&version=1.1.0&typename=PUBLIC:LBGEBRPERC2019&srsname=EPSG:3857&outputFormat=application/json&count=1000&bbox=" +
      extent.join(",") +
      ",EPSG:3857";

This is the link to my cordessandbox project.
https://codesandbox.io/embed/vector-wfs-7thkt?fontsize=14&hidenavigation=1&theme=dark

If I change the webservice URL the problem is gone and everything behaves normaly. But the point is to get it done with this URL.
If you are experiencing problems with CORS using the standard URL, you can install a CORS extension in Chrome from this link: https://chrome.google.com/webstore/detail/moesif-orign-cors-changer/digfbfaphojjndkpccljibejjbppifbc?hl=nl

var text =
      "https://geoservices.landbouwvlaanderen.be/PUBLIC/wfs?service=WFS&request=GetFeature&version=1.1.0&typename=PUBLIC:LBGEBRPERC2019&srsname=EPSG:3857&outputFormat=application/json&count=1000&bbox=" +
      extent.join(",") +
      ",EPSG:3857";

import "ol/ol.css";
import Map from "ol/Map";
import View from "ol/View";
import GeoJSON from "ol/format/GeoJSON";
import { Tile as TileLayer, Vector as VectorLayer } from "ol/layer";
import { bbox as bboxStrategy } from "ol/loadingstrategy";
import BingMaps from "ol/source/BingMaps";
import VectorSource from "ol/source/Vector";
import { Style, Fill, Stroke } from "ol/style";
import Select from "ol/interaction/Select";

var raster = new TileLayer({
  source: new BingMaps({
    imagerySet: "Aerial",
    key: "AsVf4lj-tiANM5N4_P56DC_oQQM9fjb0lMosBxFtgovzGEgcMnQuqYpeKpX-1KL2"
  })
});

var vectorSource = new VectorSource({
  format: new GeoJSON(),
  url: function (extent) {
    var text =
      "http://localhost:3030/map/https://geoservices.landbouwvlaanderen.be/PUBLIC/wfs?service=WFS&request=GetFeature&version=1.1.0&typename=PUBLIC:LBGEBRPERC2019&srsname=EPSG:3857&outputFormat=application/json&count=1000&bbox=" +
      extent.join(",") +
      ",EPSG:3857";

    console.log(text);
    return text;
  },
  strategy: bboxStrategy
});

var vector = new VectorLayer({
  source: vectorSource,
  style: new Style({
    fill: new Fill({
      color: "rgba(255,0,0,0.3)"
    }),
    stroke: new Stroke({
      color: "rgba(255,0,0,0.3)"
    })
  })
});

var map = new Map({
  layers: [raster, vector],
  target: document.getElementById("map"),
  view: new View({
    center: [589973.4805179046, 6575521.818939352],
    maxZoom: 19,
    zoom: 15
  })
});
// ref to currently selected interaction

var changeInteraction = function () {
  var select = new Select();
  map.removeInteraction(select);
  map.addInteraction(select);
};

/**
 * onchange callback on the select element.
 */
changeInteraction();


texts: Getting max value from field in Field Calculator of ArcMap

texts:
I'm working from a geodatabase in ArcMap 10.5, and using Model Builder.

I currently have a field with values that range from 0 to 30 and I want to convert that to a 0-10 scale. The idea was to use a formula like:

([VALUE]/[MAX VALUE])*10

Which, after some minor research, I assumed would be something like this in Field Calculator (using...):

(INDEX/max(INDEX))*10

But it doesn't work. The calculator does not return any error, it just doesn't give the exepected results. I tried adding [] to the max function, adding some extra fields so i has something else to search on, etc, but I always get one of these results:

All with value = 1 and some NULLs (which is odd 'casue there are no initial NULL values)
All with value = 10 and some NULLs (same as above in this case)

Any of these makes sense taking the formula in mind, but there is also no error message, so I assume there is some piece I am missing to make the calculation work correctly.

How can I (if it's even possible) get the max value in a field to calculate this formula?

Doesn't matter if it's VB Script or Python (currently using the latter) or even if it's some workaround, as long as it works.


texts: Can I just swap out a .prj file to reproject a shapefile?

texts:
I'm an experienced software developer -- but a GIS noob.

I have a shapefile that I am trying to turn into a SVG map with Kartograph. When I try to run Kartograph I get this error:

ERROR 6: No translation for Lambert_Conformal_Conic to PROJ.4 format is known.

Because I don't want to delve into Kartograph, it seems like I need to reproject the shapefile. The Wikipedia page for shapefiles says that the .prj file describes the projection.

.prj — projection format; the coordinate system and projection information, a plain text file describing the projection using well-known text format

Can I just swap out the .prj file to reproject the file?

If so, where can I find a .prj file that will work to reproject the map so that Kartograph can turn it into an SVG?


texts: Loop through coordinates, adding a marker to the map for each using HERE maps

texts:
I am trying to add a number of markers to my map, using HERE maps API. I have my coordinates in a list, and there are 30+ pairs.
var coords = [[60.1697768639186, 24.82923586371288],
[60.170497332447866, 24.828592133548653 ],
[60.17093227431904, 24.827707004572943 ], ..........

https://developer.here.com/documentation/maps/3.1.22.1/dev_guide/topics/marker-objects.html
The documentation outlines that this is the way to add single markers
var parisMarker = new H.map.Marker({lat:48.8567, lng:2.3508});
map.addObject(parisMarker);

var romeMarker = new H.map.Marker({lat:41.9, lng: 12.5});
map.addObject(romeMarker);

var berlinMarker = new H.map.Marker({lat:52.5166, lng:13.3833});
map.addObject(berlinMarker);

So I ultimately want to loop through each element of my coordinate list and add a marker for each pair.
I am trying this..
// chargepoints = 2 if the chargepoints flask route is completed
if (chargepoints = 2) {

  // ADD THE CHARGEPOINTS TO THE MAP
  for (var i=0; i<coords.length; i++) {

    var chargemarker = new H.map.Marker({lat:coords[i][1], lng: coords[i][0]});map.addObject(chargemarker);

    map.addObject(chargemarker);

}
}

However the markers are not appearing. When I insert some consolg.log("hello"); into different stages in the loop, they appear in the console however.


texts: Editing Cursor for GPS live tracking in QGIS 2.4

texts: I am interested in using QGIS 2.4 for vessel navigation due to its ability to display NOAA vector charts or BSB's as well as user created shapefiles.  I currently oversee a project that requires seafloor object buffers to be displayed so that the vessel can navigate around them during bottom trawl surveys to avoid snagging the net.  I would like to be able to modify the displayed cursor to include a steer whisker akin to what you would see in a vessel chart plotter.  Does anyone have any suggestions on how to go about this?


texts: Default folder structure in ArcGIS for Server?

texts:
I anticipate having to develop a front-end for an ArcGIS server, and until it's ready I've just been poking around the web interface of an existing instance to get a general understanding of where things are.

For example, I can find the main list of available services here, and the generic API documentation here.

What other standard locations exist in a default ArcGIS setup, if no customization is done (or, equivalently, where I can find such information)?

I looked into the links posted by @Mapperz in the comments, and while they were informative about migrating to 10.2, they didn't talk much about conventions or defaults for the ArcGIS web interface. The second link mentioned the arcgis/rest/services convention, but that was one of the ones I had already discovered.


texts: UTM coordinate conversion to Geographic coordinates

texts:
I have about 500 separate files containing GPS points in UTM coordinates from multiple zones.

I am looking for a simple and fast python or Php script that i can use to carry out a batch conversion for these files.

Thanks


texts: Field values round off while exporting XY data to a shapefile

texts:
I am trying to export a excel sheet to a shapefile where I am performing the below shown steps:

Step 1: I right clicked on the excel sheet which I wanted to export into a shapefile and selected Display XY Data as seen in the image

Step 2: After this step a Sheet$Events layer was added in my TOC.
On which I right clicked and performed export data to export it to a shapefile.

Step 3: After opening the attribute table, I have a mismatch in the values as they are rounding off.

I don't want these values to round off. It is affecting my analysis. I can't find a way to solve this problem. Please help!


texts: Model to create layer from joined MS access table

texts:
I am using arc10. I created some queries in MS access that I can open up in arcmap utilizing the database connections. It works fairly well but is rather slow.
I wanted to know if there is a model that I could create, such that when arcmap is started, it would take the joined table from access and create a layer for the user to use as I think the connection from MS access query to arcmap 10 is making the process very slow. It works but if the table were saved as a layer it is faster. I could manually do it but if there is a way to automate the process that would be better as multiple users are involved.
I'll try creating the model and seeing if I can get it to work.


texts: Converting Google Maps API to DOJO/ArcGIS API for JavaScript?

texts:
I have a GIS webpage in HTML5/JavaScript using the Google Maps API. I found out we can no longer use the Google Maps imagery provided, and I need to switch to ESRI Imagery (free).

Can someone help me is how I go about this?

As I understand the ESRI stuff uses DOJO and since I am already using JavaScript can I simply change some of my code to work?

I have a map with a few markers, lines and a circle polygon (around the marker) nothing too crazy.


texts: How to display drawings from ArcGIS on Google Maps using JavaScript or PHP?

texts:
my client has ArcGIS installed on his PC
 and he has designed roads, hospitals, houses, and parks of a specific area in ArcGIS
 he needs his system online.

challenging task for me is how can I get all drawing information from his ArcGIS that I can show that same drawing on google map using javascript or php or whatever
 I searched it but I am not satisfied yet.

there is also a solution of creating ArcGIS Server, but that will be my 2nd priority.
 I am thinking to somehow export all drawing information from ArcGIS and feed it to google map.
 bcz there will be no frequent changes in drawing he has. (means drawings are static and not changing frequently.)


texts: How do I clip buffers to multiple lines/polygons?

texts:
I have a map with 6000+ 1km buffers of ponds I want to clip them so the buffer does not cross a barrier which is a series of polylines in a separate shapefile.

I have created a 1 km buffer around all ponds.

I have merged several road files to create a barrier file.

I now want to clip the buffer rings against the buffers.

I tried converting the barrier lines to polygons and then clipping but this just clipped the edge of the whole shape not to the individual road features.

I basically do not want the buffers to cross a barrier.


texts: Using geometry generator to create line hatch fill

texts:
Since QGIS rasterizes any sort of pattern fill, I have been searching for a way to create, through the geometry generator, a line-hatch fill for a polygon. I so far haven't found anything I can easily adapt to the purpose. I found this for a point fill, which is similar to what I'd like for a line-hatch fill: https://gis.stackexchange.com/a/337774/122005. I don't care about random spacing, at least at this point. All I care about are two parameters: line spacing and angle. Everything I have found is much more involved and I can't pare the code down to just what I need.
Here's an image of what I'm after:

Many thanks in advanced for anyone's assistance.


texts:
I'm working on a script to replace a layer datasource, and I'm getting an error I cannot figure out.

When I run the following code in a python window inside ArcMap, I get a "Runtime error : 'unicode' object is not callable"

import arcpy
mxd = arcpy.mapping.MapDocument ("CURRENT")
for lyr in arcpy.mapping.ListLayers(mxd):
     if lyr.supports("DATASOURCE"):
         if lyr.dataSource (r"N:\3_PLANNING\6_GIS\GIS Data\INDIANA GIS\Indiana_County.mdb\Boundaries\Indianacounty_shp"):
             lyr.replaceDataSource (r"N:\3_PLANNING\6_GIS\GIS Data\Indiana_County_Database_201201.gdb", "FILEGDB_WORKSPACE", r"N:\3_PLANNING\6_GIS\GIS Data\Indiana_County_Database_201201.gdb\Boundaries\County_Boundary")

Problem is, I'm not very good at python yet, and I'm not sure what the unicode object is that the error is referring to.


texts: Query slowed down by gist-index on points

texts:
So I have this simple query, selecting all points in my table that are within a given polygon:

WITH my_polygon as (
  SELECT ST_GeomFromEWKT('<EWKT polygon>') as geom
)
SELECT count(*)
FROM points as a
,    my_polygon as p
WHERE ST_Within(points.geom, my_polygon.geom);

I have two equal points tables, both consisting of the same kind of points, scattered globally. The only difference is:

One is sparse _____ 500.000 rows
One is dense __ 340.000.000 rows

I have now executed the query on both tables, testing with and without gist-index on the points.geom column. On the sparse table I get these execution times:

No indexing ________ 483.544 ms
With gist index ________ 0.142 ms

Looks good, index is working properly. Now for the dense table:

No indexing _____ 195307.138 ms
With gist index ___ 234495.684 ms

This was not as expected. Is the gist index not supposed to speed up this query?
The explain analyze for all cases is provided in this GitHub-Gist

(I did performed vacuum analyze after creating the index.
Also I see the same behavior when gin-indexing an array. Query slows down on the big table using the index.)


texts: How to compute (X,Y,Z) of a point from given attributes

texts:
Is there a way to compute the geodetic coordinates (Earth Centered) of a satellite S(x,y,z). The original problem is to compute the footprint of an antenna,
given the following information:

Antenna(x,y,z) geodetic coordinates (Earth Centered) for an antenna
e = antenna elevation [0:90]
Az = (antenna azimuth) [0:360]
R = Radius of Earth (appr. 6378.1363 km)
h = height of S (in kilometers)

If it turns the solution easier, the coordinates of A and S can be written in (lat,long,h) format:

Antenna(lat,long, h = 0)
S(lat,,long,h)

I'm trying to write the code in python, so I can use some libs (like pylint and pyproj, for instance)


texts: Angle orientation for data defined label rotation in QGIS?

texts:
It seems that angle orientation for data defined label rotation is different in QGIS 2 (anti-clockwise) and QGIS 3 (clockwise). I mean angle values that are stored in the defined field when modifying label orientation with the graphic tool of the label toolbar.
Is this intentional ?


texts: Converting a layer file (*.lyr) to a QGIS compatible file without using ArcGIS

texts:
I have some .lyr files and I need to be able to open them in QGIS.
I don't have access to ArcGIS so I need to convert those files to another format compatible with QGIS (.qlr, .xml, .sld, .style,...) either using QGIS or another tool (preferably an open source one) .


texts: Magellan GPS stops recording data in suspend mode

texts:
Currently having a problem with my Magellan Explorist 510 GPS unit. After I start recording the track on the unit, after some time, it just stops recording data in suspend mode (which should turn everything off except the tracking of GPS position on the track).

I don't know if it's a setting that I have changed or if it's a problem with the unit. Before I try a firmware reset, I was just wondering if anyone else has experienced this issue?


texts: select feature info on user defined polygon

texts:
Ok so I have a Openlayers App which is displaying few of the features and I now want to give the user the freedom to create a Box or polygon. All the features which are overlapped by this user defined box, I wish to export the data and then display it onto a table in the html.

Is there any OL example to such application requirement?
Which would be better to extract the feature info, WFS or WMS?

I know it can be done in GMaps but I have OSM and GMaps both as base maps so I require a OL functionality.


texts: Extend / override Mapbox-gl 'baselayer' loadTile

texts:
Similar to how Leaflet's guide shows how to override L.TileLayer's getTileUrl method in place of your own (and L.GridLayer's createTile method),

I want to know how to extend/override Mapbox-gl's baselayer tiles (as opposed to only using a tile server URL to add in custom tiles).

I'm looking to avoid just adding a layer and instead want to overwrite the baselayer itself.

I saw this issue discussion of customizing RasterTileSource.loadTile and wondered if this was the right direction.

Leaflet's example:

L.TileLayer.Kitten = L.TileLayer.extend({
    getTileUrl: function(coords) {
        var i = Math.ceil( Math.random() * 4 );
        return "https://placekitten.com/256/256?image=" + i;
    }
});

L.tileLayer.kitten = function() {
    return new L.TileLayer.Kitten();
}

L.tileLayer.kitten().addTo(map);


texts: Multiple ring buffer (MRB) with different buffer values according to 2 fields

texts:
I have a shapefile with farms information.
I know that the size of this farms is increased during last year by a certain percentage. For each row I calculate (according to farm area) the buffers radius pre- and post-growth.

So now I have two fields with the buffer radius values (that are different  for each row).

How can I run a MRB according to these fields value?


texts: ArcGIS Online fails to add zipped File Geodatabase?

texts:
I couldn’t figure out why The ArcGIS Online fails to add a zipped “file geodatabase”, and gives the following message:

  Invalid name ‘xxx.gdb’ found in input publish parameters


texts:
I have an Enterprise geodatabase and want to connect to the MSQL database from QGIS. I have been trying to connect to MSQL, but the connection couldn't be done as I am getting an error message. I am getting the following error:

Error opening connection: [Microsoft][ODBC SQL Server Driver]Cannot generate SSPI context [Microsoft][ODBC Driver Manager] Driver's SQLSetConnectAttr failed QODBC3: Unable to connect


texts: How to open QGIS project properties focused on "CRS" and select choosen projection

texts:
I'm currently working on QGIS plugin as part of a school project.  One feature of this plugin is going to be the possibility to transform coordinates from one EPSG code to another.

Therefore it would be the most useful and user-friendly way to let the user open the "project properties" dialog by clicking a specific button in the plugin, in this dialog he chooses the CRS he wants to use and on closing the dialog the chosen value is taken to the plugin code so I can use it for the transformation itself.

I already searched in the QGIS API but without any success. I made it to open some other menus via the "iface" reference in my code, but unfortunately not the one you open by clicking in QGIS at to bottom right corner, left of the bubble which opens the message view.

Hope I made my problem as understandable as possible, if not feel free to ask more details.


texts: Seeking GIS data on military activity in UK, i.e. low flying zones, military practice areas, etc?

texts:
I'm trying to source shapefiles of low flying zones, military aerodromes, military practice areas for the Wales area including offshore in the 12nm limit and EEZ.

Where could I source this information?


texts: Creating shapes/lines from geo coordinates in postgis

texts:
Following the suggestion from this topic I'm posting my question here.

I have a table in postgres/postgis with geo data and I want to show them on map via qgis.

Here is the table:

city  |  street  |  geo
A     |  A1      |  LINESTRING(17.257125 54.595952,17.253573 54.599437)

I do not know how to do this. I can show a single point on map (creating geo column with AddGeometryColumn function) without any problem, but how to show this part LINESTRING(17.257125 54.595952,17.253573 54.599437) (points with line between them). Any advices, suggestions or link to some articles about it?

I was thinking about splitting geo coordinates from Linestring(), but in the end this will not work at least that is my opinion.

I'm new to GIS things so please guide me.


texts: Add svg filter to paths created by leaflet?

texts:
If I were to create a simple inline svg, I could modify paths to have a drop shadow or a blur by creating a <filter> within a <defs> grouping. However, Leaflet dynamically generates the SVGs it uses to display paths, and I'm not sure how to add a filter after the fact.

Failed Ideas:

I've tried using jQuery to append a defs to the SVG after the fact:

Like this:

    $("svg").prepend("<defs><filter id='dropshadow' height='130%'><feGaussianBlur in='SourceAlpha' stdDeviation='3'/> <!-- stdDeviation is how much to blur --> <feOffset dx='2' dy='2' result='offsetblur'/> <!-- how much to offset --><feMerge><feMergeNode/> <!-- this contains the offset blurred image --><feMergeNode in='SourceGraphic'/> <!-- this contains the element that the filter is applied to --></feMerge></filter></defs>");

And then specify the location of the defs in my CSS, like this:

    filter: url("../index.html#dropshadow");

Have the <defs> / filter definition already inline in the html above the map element, and then reference the filter definition from the svg via a the CSS filter URL.

My code (inline above the leaflet map div):

    <svg xmlns="w3.org/2000/svg"; version="1.1">
        <defs>
            <filter id='dropshadow' height='130%'>
                <feGaussianBlur in='SourceAlpha' stdDeviation='3'/> <!-- stdDeviation is how much to blur -->
                <feOffset dx='2' dy='2' result='offsetblur'/> <!-- how much to offset -->
                <feMerge><feMergeNode/> <!-- this contains the offset blurred image -->
                <feMergeNode in='SourceGraphic'/> <!-- this contains the element that the filter is applied to --></feMerge>
            </filter>
        </defs>
    </svg>


texts: Calculating water occurrence of Sentinel-1 images in Google Earth Engine

texts:
I would like to derive the occurrence of surface water bodies from Sentinel-1 image collection in the Google Earth Engine.In my example I could already compute the water mask using Otsu’s method for image segmentation: https://code.earthengine.google.com/499580510635aa5f421886d4e409af03.

This procedure has now to be applied to the entire image collection, so that the variable threshold is recalculated for each image and then the water mask is determined.The occurrence of water bodies should be presented relatively (0%-no occurrence, 100%-always water present) as shown in Global Surface Water Explorer.

Does anyone have a solution or a hint for this?

Create the S1 image collection:

var s1 = ee.ImageCollection('COPERNICUS/S1_GRD')
  .filterBounds(Map.getBounds(true))
  .filterDate('2019-01-01','2019-12-31')
  .filter(ee.Filter.listContains('transmitterReceiverPolarisation', 'VV'))
  .map(function(image){return image.clip(Map.getBounds(true))});
print(s1);

Otsu's method for image segmentation:

var histogram = image.select(s1_band).reduceRegion({
  reducer: ee.Reducer.histogram()
      .combine('mean', null, true)
      .combine('variance', null, true),
  geometry: polygon,
  scale: 10,
  bestEffort: true
});
print(histogram);

// Chart the histogram
print(Chart.image.histogram(image.select(s1_band), polygon, 10));

// Return the DN that maximizes interclass variance in S1-band (in the region).
var otsu = function(histogram) {
  var counts = ee.Array(ee.Dictionary(histogram).get('histogram'));
  var means = ee.Array(ee.Dictionary(histogram).get('bucketMeans'));
  var size = means.length().get([0]);
  var total = counts.reduce(ee.Reducer.sum(), [0]).get([0]);
  var sum = means.multiply(counts).reduce(ee.Reducer.sum(), [0]).get([0]);
  var mean = sum.divide(total);

  var indices = ee.List.sequence(1, size);

  // Compute between sum of squares, where each mean partitions the data.
  var bss = indices.map(function(i) {
    var aCounts = counts.slice(0, 0, i);
    var aCount = aCounts.reduce(ee.Reducer.sum(), [0]).get([0]);
    var aMeans = means.slice(0, 0, i);
    var aMean = aMeans.multiply(aCounts)
        .reduce(ee.Reducer.sum(), [0]).get([0])
        .divide(aCount);
    var bCount = total.subtract(aCount);
    var bMean = sum.subtract(aCount.multiply(aMean)).divide(bCount);
    return aCount.multiply(aMean.subtract(mean).pow(2)).add(
           bCount.multiply(bMean.subtract(mean).pow(2)));
  });

  print(ui.Chart.array.values(ee.Array(bss), 0, means));

  // Return the mean value corresponding to the maximum BSS.
  return means.sort(bss).get([-1]);
};

var threshold = otsu(histogram.get(s1_band+'_histogram'));
print('threshold', threshold);

Note on a tripping hazard for calculating the occurence:

The .filterBounds method based of an area of interest also selects truncated satellite images. Unfortunately I have no solution to prevent this.

Example of a truncated image


texts: How to remove layer from layer control in Leaflet

texts:
I have an AJAX call to dynamically create heatmap from database on an addOverlay layercontrol. This event is handled by a button whose text changes upon click event. I am able to add layer to map but not able to remove (when I remove and add again by clicking "generate heatmap", a new "HWC HeatMap" is added each time on LayerControl). How can I solve this problem?
<button id="btnHeatMap" class='btn btn-primary btn-block'>Generate HeatMap</button>

$("#btnHeatMap").click(function(){
    var heatMapPoints = [];
    var heatHWC;
    $.ajax({
        url:'heatmap.php',
        success: function(response){
            jsnHWC = JSON.parse(response);
            jsnHWC.features.forEach(function(feature){
                heatMapPoints.push([feature.geometry.coordinates[1],feature.geometry.coordinates[0], 0.5]);
            });
            heatHWC = L.heatLayer(heatMapPoints, {
                radius: 25,
                minOpacity: 0.4,
                gradient: {0.3: 'green', 0.7: 'yellow', 1: 'red'}
            });

            if($("#btnHeatMap").html()=="Generate HeatMap") {
                heatHWC.addTo(mymap);
                ctlLayers.addOverlay(heatHWC, "HWC HeatMap");
                $("#btnHeatMap").html("Remove HeatMap");
            } else {
                if(mymap.hasLayer(heatHWC)){
                    mymap.removeLayer(heatHWC);
                }
                $("#btnHeatMap").html("Generate HeatMap");
            }
        },
        error: function(xhr, status, error) {
            console.log("ERROR: "+error);
        }
    });
});


texts:
I have a problem to rotate a point (valve symbol) on a polyline I used the code below and it works with me but when I put point on line only and did not work with polyline

 public double pointAngle(IPoint aPoint, IFeature LineFeature)
    {
        IPolyline pPolyline = LineFeature.Shape as IPolyline;
        IPoint pFromPoint = pPolyline.FromPoint;
        IPoint pToPoint = pPolyline.ToPoint;

        ILine pLine = new Line();
        pLine.PutCoords(pToPoint, pFromPoint);
        ILine pLine2 = new Line();
        pLine2.PutCoords(aPoint, pFromPoint);
        ILine pLine3 = new Line();
        pLine3.PutCoords(aPoint, pToPoint);
        double dAngle;
        if (pLine3.Length < pLine2.Length)
        {
          //  dAngle = ((pLine3.Angle* 180)/pi); //Math.PI) ;
            dAngle = ((pLine3.Angle * 180) / Math.PI) ;
        }
        else
        {
            dAngle = ((pLine2.Angle* 180)/ Math.PI) ;
        }
        dAngle = dAngle + 90;
        return dAngle;

    }

and when I used this code nothing happen with me because I think I use ITransform2D.rotate wrong

 IFeature feature = dssClass.SelectFeatureByID("Valve", 180);
        IPoint ppPoint = feature.Shape as IPoint;
        double radian = 1.5708;
        ITransform2D pTrans2D = ppPoint as ITransform2D;
        pTrans2D.Rotate(ppPoint, radian);

All I need like this picture


texts:
When I try to use the following expression, the preview shows me the correct sum, but when I click ok I get this notification. This tells me it can't calculate column "Subsidie J"

It seems that it can't find the value for column "Subsidie J" but it's all in the attribute table.

Am I using the wrong expression? The value of "Subsidie J" is a CASE IF .... THEN ... END expression. Now I want the sum of these values sorted by owner.


texts: Error in passing geometry column to SQL query

texts:
When I execute the following query with a static table name it's working fine:

SELECT pe.structure_category_name,SDO_NN_DISTANCE(2),pe.bt_object_subid,
pe.bt_owning_node_id INTO v_sDPsct,v_saddDPdist,v_dpTypeNearSrvAdd,
v_sbtOwningNodeID FROM piper_equipment pe WHERE
SDO_NN(pe.shape,(unservedValue.shape),'sdo_batch_size=10',2) = 'TRUE' AND ROWNUM < 2 AND pe.bt_object_id = v_saddDPid

Here unservedValue is the data_record fetching through cursor.
I have to pass table name as variable , so I have changed above query into the following:

stmt_piperequipment3 := 'SELECT pe.structure_category_name,SDO_NN_DISTANCE(2),
pe.bt_object_subid,
pe.bt_owning_node_id INTO v_sDPsct,v_saddDPdist,v_dpTypeNearSrvAdd, v_sbtOwningNodeID
FROM '||VWPIPERDATA||' pe WHERE
SDO_NN(pe.shape,(unservedValue.shape),''sdo_batch_size=10'',2) = ''TRUE'' AND
ROWNUM < 2 AND pe.bt_object_id = v_saddDPid';
execute immediate stmt_piperequipment3;

While executing above query it's giving following error

ORA-00904: "UNSERVEDVALUE"."SHAPE":invalid identifier.

My question is how to pass a shape into the string variable?


texts: My data at a certain latitude is disappearing in Mapbox Studio?

texts:
I'm trying to make a simple map showing 4 villages in Alaska (and Fairbanks for reference) but for some reason the northernmost village is being cut-off it seems. When viewing the data in Mapbox Studio there is a mysterious line that seems to be the culprit? I have no idea what this line is from. See picture below.

The coordinates for the disappearing village are -152.65091, 66.56305. Here is the data file in geojson format as a gist villages.geojson gist

And here is the working map Map


texts: Problem finding the right CRS for CSV file

texts:
I got a problem regarding a CSV layer I have to use. I've imported and tried different import setting still everytime QGIS shows the points at the African coast instead of Germany. The CSV has two attributes for lat/lng but I can't tell if there might be something wrong with them.

For visualization here the screenshot of the data and import setting:

And here is the outcome after importing to the map:

Can someone give some advice on what might be wrong?


texts: Transforming geostationary satellite image to lon/lat

texts:
The ICARE Data and Services Center distributes the SEV_AERUS-AEROSOL-D3 product in HDF5 format (http://www.icare.univ-lille1.fr/archive/?dir=GEO/MSG+0000/SEV_AERUS-AEROSOL-D3.v1.03/).

The main component of the HDF file is a 3712 x 3712 matrix, which represents measurements on a geostationary view of the earth, such as in the following image:

In R, I am able to read the latter into a matrix object named m using the h5 package:

library(h5)
setwd("/home/michael/Dropbox/BGU/Alex/HDF5")
filename = "SEV_AERUS-AEROSOL-D3_2005-01-01_V1-03.h5"
h = h5file(filename)
m = h["ANGS_06_16"][]

Here are the matrix dimensions:

dim(m)
## [1] 3712 3712

And, for example, these are the values in the first 1-5 rows/columns:

m[1:5, 1:5]
##        [,1]   [,2]   [,3]   [,4]   [,5]
## [1,] -32768 -32768 -32768 -32768 -32768
## [2,] -32768 -32768 -32768 -32768 -32768
## [3,] -32768 -32768 -32768 -32768 -32768
## [4,] -32768 -32768 -32768 -32768 -32768
## [5,] -32768 -32768 -32768 -32768 -32768

The HDF5 file also provides metadata with the lon/lat of two corners for the bounding box of m (below) as well as the projection type (specified simply as "GEOS"). Here is a reproduction of the bounding box coordinates exactly as given in the metadata:

library(raster)
ext = new("Extent"
          , xmin = -75.1001529854046
          , xmax = 75.1001529854046
          , ymin = -77.9875583358503
          , ymax = 77.9875583358503
)
ext
## class       : Extent
## xmin        : -75.10015
## xmax        : 75.10015
## ymin        : -77.98756
## ymax        : 77.98756

The problem is that I need to find out the spatial location of each measurement in the matrix. In order to do that, the matrix needs to be converted into a point layer with lon/lat coordinates for each pixel.

My idea was:

Convert the latter bounding box coordinates of m to "GEOS" projection coordinates.
Convert m to a RasterLayer object, using the "GEOS" projection bounding box and a "GEOS" CRS definition.
Transform the resulting RasterLayer to a point layer.
Reproject the point layer to lon/lat.

Unfortunately, I am stuck in step (1). For example, the following code takes the Southernmonst/Westernmost corner in order to calculate its coordinates in the "GEOS" system:

ext = data.frame(x = ext[1], y = ext[3])
ext
##           x         y
## 1 -75.10015 -77.98756
coordinates(ext) = ~ x + y
proj4string(ext) = '+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0'
spTransform(ext, "+proj=geos +lon_0=0 +h=35785831 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs")

The last expression gives the following error message:

Error in spTransform(x, CRS(CRSobj), ...) :
error in pj_transform: tolerance condition error

In other words, the (-75.10015, -77.98756) point can't be converted from lon/lat to "GEOS". I also tried different variations of the "GEOS" PROJ.4 string, none worked. The "tolerance condition error" is apparently a general PROJ.4 issue, since I also get the same error when trying to reproject the (-75.10015, -77.98756) point from lon/lat to "GEOS" using pyproj in Python, as suggested here.

When using coordinates closer to the equator, such as (10,20), the code works fine (no "tolerance condition error"). But this doesn't help, since I can't subset the m matrix without knowing the new bounding box coordinates.

Another option could have been translating the matrix to lon/lat using Longitute/Latitude static files, as suggested here. However I couldn't find any such files matching the 3712 x 3712 dimensions.

Does anyone have experience with going from a 3712 x 3712 matrix representing a geostationary satellite image to a point layer in lon/lat? Are there any static files providing predetermined lon/lat values of each cell in such a grid?


texts:
In attempting to find at least a partial solution for Google Earth font change, changing the font for my own placemark balloons, I ran into another problem.

If I use <BalloonStyle> to give styling for balloons then every placemark affected by this style is given a balloon whether or not there is any $[description] text.  Normally the balloon is suppressed if the description is empty.

For example I may wish to use:

<BalloonStyle>
  <text>
    $[description]
  </text>
</BalloonStyle>

However if I use this style for a placemark without a description I get an empty balloon.

Can I give styling using <BalloonStyle> and still have suppression of balloons for empty descriptions?


texts:
I try to download data from http://download.geofabrik.de/europe for UK. I downloaded England, Wales and Scotland shapefiles. I have also Ireland and Northern Ireland shapefile (http://download.geofabrik.de/europe/ireland-and-northern-ireland.html) but I don't know how to cut them to obtain only Northern Ireland.

Attributes of all shapefiles (England, Wales and Scotland, Ireland and Northern Ireland) look like:

.

Can someone help me?


texts: Error with ST_Value, ST_NearestValue: out of range raster coordinates

texts:
I have a UK 5m DTM coverage loaded into Postgis with:

raster2pgsql -s 27700 -I -C *.tif -t 100x100 ....

which looks more or less like the below for a full 1000 x 1000 pixel tile, before tiling:

This works as expected, however when I run a query for ST_Value, such as:

SELECT ST_Value(rast, geom, true)
  FROM terrain_5, addresses
 WHERE ST_Intersects(rast, geom)
   AND address_id = xxxx;

I get the error:

  NOTICE:  Attempting to get pixel value with out of range raster coordinates: (2, 100)

If I run the same query, but using ST_NearestValue instead I get two rows back, eg,

-3.40282346638529e+38
 70.7099990844727

I have checked for the points where this happens and they are always round numbers divisible by 500.

If instead I request the metadata, ie, the same query but using (ST_MetaData(rast)).* I get two rows,

upperleftx | upperlefty | width | height | scalex | scaley | skewx | skewy | srid  | numbands

166500 |      27500 |   100 |    100 |      5 |     -5 |     0 |     0 | 27700 |        1

166500 |      27000 |   100 |    100 |      5 |     -5 |     0 |     0 | 27700 |        1

It is pretty obvious that the problem is due to 100x100 tiles of pixel size 5m, resulting in edge conditions every time a coordinate being tested for intersection is divisible by 500.

I have solved the problem by rewriting the above query using max(ST_NearestValue(rast, geom)), ie,

SELECT max(ST_NearestValue(rast, geom, true)
  FROM terrain_5, addresses
 WHERE ST_Intersects(rast, geom)
   AND address_id = xxxx;

the final parameter makes no difference, as there are no no data values in the coverage, as can be seen from the input raster.

If I run the histogram, excluding, no data, I get and add up the values:

SELECT sum(count)
  FROM
      (SELECT
            (ST_Histogram(rast, 1, true, 10)).count
         FROM terrain_5, addresses
        WHERE ST_Intersects(rast, geom)
         AND address_id = xxxx
      ) foo;

I get 20000, ie, the sum of the pixels in both the intersected tiles that intersect the point divisible by 500.

I have seen this question, but I don't think it really answers my case.

I have a solution, using max(ST_NearestValue(rast, geom)), but this does strike me as unexpected behaviour, or I am missing something really obvious. Sadly, when I ran this on the national dataset, the value that was chosen from the two results was the low one, so I have a lot of houses whose height above sea level is -3.40282346638529e+38, when even by British standards of building in flood plans, strikes me as excesssive. I only noticed when looking at summary statistics after the fact, as my original update query swallowed the errors.


texts:
We have an animation built with leaflet that was working until about a month ago.

Now, the play button is just like an image, and nothing happens when I hover or click on it.
In the inspect element under sources I have the following:

Could not load content for https://cdnjs.cloudflare.com/ajax/libs/vis/3.12.0/dist/vis.js (HTTP error: status code 404, net::ERR_HTTP_RESPONSE_CODE_FAILURE)

May this be the cause and if yes, how I may fix it?


texts:
I have QGIS 2.6.1 and when I want to get some results using gis2web (by Tom Chadwin) plugin, in fact, to get a webmap in OpenLayers, it seems that plugin doesn't work properly.
E.g. I have a point layer and when I exported with this plugin, the points' labels are overlapped on them.
Another thing that is "Highlight features" and "Show popups on hover" don't work either.
Had anyone any kind of this trouble with this plugin?
PS. I have gis2web 0.21.3 version.


texts: Postgres Layer not showing as editable when creating an offline project

texts:
My postgres layer in a QGIS saved project does not appear in the "Creat Offline Project" dialog as editable.  I want to check it out for offline editing?  If I uncheck the "Show only editable layers" box, it then appears.

It checks out fine but edits cannot be saved in the remote session.


texts:
I have a direct broadcast AQUA data. I used RT-STPS to extract the level 0 .PDS files and process them further to level 1A and level 1B using the modisL1DB software. These produce .hdf files.

I can view the images however, they are not geo-referenced. How can I geo-referenced these data?


texts: postgres / postgis detecting new roads algorithm

texts:
I have a postgres (psql 9.1.22) database containing all the lines that ever existed in a region in OSM. I would like to detect newly mapped roads by year. This has to be done by a spatial query (because of how OSM works). The algorithm below is based on previous similar analysis in QGIS and FME. It works in theory. But I would like to hear if this way of approaching things is efficient in Postgres.

The idea is to get a collection of just the new roads for every year and use it for some stats (dump) and visuals (QGIS).

Reproject to a local projection, so distances are in meters
Extract valid roads at 1/1/2007 (year 0 for practical purposes) to  view timeslice0
In timeslice0: calculate variable newlength = STLength (all roads were new). Copy timeslice0 to timeslicecollection
Buffer timeslice0 by 7 meters
Create timeslice1 with roads valid 1/1/2008
Intersect timeslice1 on top of timeslice0, results to timeslice1intersect
In timeslice1intersect, calculate length of pieces outside of buffers into variable newlength
Aggregate newlength by osm_id (original line id) and pass the value on to timeslice1.
Add timeslice1 to timeslicecollection
Repeat for timeslice2 over timeslice1, till timeslice for 2017
In timeslicecollection, compute newroad=1 if newlength>7 & newlength/length > 0.25
In timeslicecollection, compute improvedroad=1 if newlength>7 & newlength/length =< 0.25


texts: ArcMap Addin - add interoperability connection programatically

texts:
I'm currently building an add-in for ArcMap - ArcView using the .Net framework.

The user of ArcMap and the Add-in should access data published from WFS services.

I am aware that it is possible to manually add WFS services via the Interoperability Connections interface.

But, is it possible to add these WFS services programatically? (chances are - my users will find it difficult to add these manually)

And if it is, how could it be done?

I can't seem to find any info on this subject.


texts: StopIteration: iteration not started error in arcpy.da.UpdateCursor

texts:
I cannot figure out why I keep getting this "StopIteration: iteration not started"error. I have confirmed the feature class 'vriTHLB_eliminated' exists and has records and all required fields.

Here is a simplified version of my code:

fieldList = ["SPECIES_CD_1","area_ha"]
with arcpy.da.UpdateCursor(vriTHLB_eliminated,fieldList) as cursor:`
    for row in cursor:
        if row[0] == None or row[0] == "":  # If species 1 is null / absent / or empty string;
             cursor.deleteRow()             # delete these rows, they do not help with analysis.

        cursor.updateRow(row) #this definitely does need to indented inside the loop like this
del row, cursor


texts: Is it possible to insert note to remember the info of a project?

texts:
Sometimes I need to open an old project, and mostly I remember how it works (i.e. the structure, the formatting and so on) but sometimes not.

At other times I open a project composed by a colleague of mine and it is hard to understand how the project was conceived.

Is it possible to insert in the project a non-printable note to help to understand the project more quickly?

edit: I use QGIS 1.8 (2.0 soon) on windows 7


texts: How to compare polygons in different layers in QGIS 3.x?

texts:
To formulate my question with a concrete example: I have two multi-polygon layers of the same city.

Layer A has polygons of the city divided in parts e.g. north, south, west, east.
Layer B has smaller polygons that represent green e.g. trees, forests, parks and such.

Now I want to add a new column "city_part" in the attribute table of layer B. And likewise I'd want an either QGIS-expression- or Python-coded script, that can figure out in which polygon of layer A every single polygon of layer B lies.
So the concrete result would be a layer B in which every polygon's got a new attribute, telling whether it resides inside north, south, west or east. How can this be done?


texts: Mapstore geostory add local image from a computer

texts:
I couldn't find information wheather is it possible to add media such as image or video by uploading or by drag'n'droping?
Now it is possible only by adding href (accessible from web) as I understand.


texts: I have created the map with some point and line layers. When I print this map, some points and lines of the map looks like shifted from their original locations. The print preview doesn't show any problems. The view of the map on the screen also looks ok. Anyone knows how to solve this problem?


texts: How do I get Visual Studio to run EsriRegasm.exe as administrator?

texts:
When run Visual Studio 10 as administrator and build a project with this in my csproj:

  <Target Name="BeforeClean">
    <Exec WorkingDirectory="$(CommonProgramFiles)\ArcGIS\bin"
      Command="esriRegasm.exe &quot;$(TargetPath)&quot; /p:Desktop /u "
      Condition="Exists('$(TargetPath)')" />
  </Target>
  <Target Name="AfterBuild">
    <Exec WorkingDirectory="$(CommonProgramFiles)\ArcGIS\bin"
    Command="esriRegasm.exe &quot;$(TargetPath)&quot; /p:Desktop " />
  </Target>

It is able to do the clean (unregister) without any problems.  However it complains saying ...

  Registration failed.  Could not write to disk.

... when it does the AfterBuild.  (Notice I've removed the /s option on the command line.)

I have a bat file that can call EsriRegasm.exe without any problems when run as administrator.  However, when I run it as a non-admin I get the same error, which is expected.  Batch files run as admin results in spawned processes that are also run as admin.  Can Visual Studio do this too?

Is there some option I can add to the csproj file to make it run EsriRegasm.exe as admin?

I've already set my UAC to "never notify", as suggested by kb 39158.  And set the file properties on EsriRegasm.exe to run as administrator.


texts:
I'm trying to build QGIS 2.8.6 on Windows 7 32 bit using MSVC Express 10. The build process completes without errors, but with a bunch of warnings. When I run the qgis.exe after installation, it only shows the splash screen for a few seconds until "Starting python" appears. Then it vanishes and nothing happens.
When building python_module_qgis__gui in MSVC Express I get the following warnings. I get similar ones when building python_module_qgis_core and python_module_gis_analysis.

Generating gui/sip_guipart0.cpp, gui/sip_guipart1.cpp,
gui/sip_guipart2.cpp, gui/sip_guipart3.cpp    sip_guipart0.cpp
c:\Program Files\Microsoft Visual Studio 10.0\VC\include\math.h(162):
warning C4211: Nonstandard extension used : redefined extern to static
sip_guipart1.cpp c:\Program Files\Microsoft Visual Studio
10.0\VC\include\math.h(162): warning C4211: Nonstandard extension used : redefined extern to static
sip_guipart2.cpp c:\Program Files\Microsoft Visual Studio
10.0\VC\include\math.h(162): warning C4211: Nonstandard extension used : redefined extern to static   sip_guipart3.cpp c:\Program
Files\Microsoft Visual Studio 10.0\VC\include\math.h(162): warning
C4211: Nonstandard extension used : redefined extern to static
Looks like there might be a problem with Python and SIP?

My cmake configuration is this:

Found PostgreSQL: C:/OSGeo4W/lib/libpq.lib
Found SpatiaLite: C:/OSGeo4W/lib/spatialite_i.lib
Found Qt version: 4.8.6
Touch support disabled
Found QScintilla2: C:/OSGeo4W/lib/qscintilla2.lib (2.8.4)
Found Python executable: C:/OSGeo4W/bin/python.exe
Found Python version: 2.7.4
Found Python library: C:/OSGeo4W/apps/Python27/libs/python27.lib
Found SIP version: 4.16.5
Found PyQt4 version: 4.11.3
Found QScintilla2 PyQt4 module: 2.8.4

Does anyone have an idea what might be wrong?


texts: Color points based on values in OpenLayers 6

texts:
I'm trying to color points based on a value called 'ML_ABSOLUTE_VALUE', in a GeoJSON dataset but I can't get it to work. Here is what I have set up so far:
Create the GeoJSON data source:
    const sourcePoints = new VectorSource({
      url: './data/dataset1_fixed.json',
      format: new GeoJSON()
    });

Next, I have a function that bins features based on the value of 'ML_ABSOLUTE_VALUE' for each feature:

      if(feature.get('ML_ABSOLUTE_VALUE') >= 0 && 'ML_ABSOLUTE_VALUE' <= 100) {
        fillcolor = '#2b83ba';
      } else if(feature.get('ML_ABSOLUTE_VALUE') > 100 && ('ML_ABSOLUTE_VALUE') <= 200) {
        fillcolor = '#91cba9';
      } else if(feature.get('ML_ABSOLUTE_VALUE') > 200 && ('ML_ABSOLUTE_VALUE') <= 300) {
        fillcolor = '#def2b4';
      } else if(feature.get('ML_ABSOLUTE_VALUE') > 300 && ('ML_ABSOLUTE_VALUE') <= 400) {
        fillcolor = '#ffdf9a';
      } else if(feature.get('ML_ABSOLUTE_VALUE') > 400 && ('ML_ABSOLUTE_VALUE') <= 500) {
        fillcolor = '#f69053';
      } else if(feature.get('ML_ABSOLUTE_VALUE') > 500) {
        fillcolor = '#d7191c';
      }

      return [new CircleStyle({
        image: new Circle({
          fill: new Fill({color: fillcolor}),
          stroke: new Stroke({color: 'black', width: 1}),
          radius: 5
        }),
      })];
    };

Last, I set up the Map() and View(), and pass in the customStyleFunction() for style.
  new Map({
    target: 'map-container',
    layers: [
    new TileLayer({
      source: new OSM()
    }),

    new VectorLayer({
      source: sourcePoints,
      style: customStyleFunction
      })],
  view: new View({
    center: fromLonLat([-70, 35]),
    zoom: 10
  })
});

I get an error 'fillcolor is not defined'. Any suggestions?


texts: QGIS - Convert a text field to a clickable hyperlink, hotlink, or URL?

texts:
Just starting out learning QGIS, on version 3.6.3 Windows 64 bit.

I have imported some data into QGIS. When I click the parcel I would like for the "WIC" field attribute to be a clickable hyperlink, hotlink, or URL. It is a link to a simple text file.

How can I do this? Here a screenshot below :


texts: Orphaned Holes after Mercator-Projection in R

texts:
In R I want to apply the Mercator-Projection to a SpatialPolygonsDataFrame before plotting it with ggplot.
I dont want to use ggplot::coord_map() since the map is supposed to be the "canvas" for a network-graph, that i dont want to be mercator-projected aswell.

After using the sp::spTransform() function I was unable to broom::tidy() the resulting SpatialPolygonsDataFrame due to "orphaned holes"...

Is the SpatialPolygonsDataFrame provided by cshapes::cshp() somehow faulty or am I doing the transformation wrong?

PS: The cleangeo::clgeo_Clean()-function (Suggested here: Fixing orphaned holes in R) has been computing for at least 3 hours sofar. It doesn't seem to be the ideal approach - even if it turns out to be working...

> library(cshapes)
Loading required package: sp
Loading required package: maptools
Checking rgeos availability: TRUE
Loading required package: plyr

> library(broom)

> wrld <-
+   cshp(as.Date("2016-01-01"))
Warning message:
use rgdal::readOGR or sf::st_read

> class(wrld)
[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

> proj4string(wrld)
[1] "+proj=longlat +ellps=WGS84"

> wrld <-
+   spTransform(wrld, CRS("+proj=merc +lon_0=0 +k=1 +x_0=0
+                          +y_0=0 +ellps=WGS84 +datum=WGS84
+                          +units=m +no_defs"))

> wrld <-
+   broom::tidy(wrld, region="COWCODE")
Error in createPolygonsComment(p) :
  rgeos_PolyCreateComment: orphaned hole, cannot find containing polygon for hole at index 84


texts: Converting LAS to DWG/DXF using open-source software?

texts:
Is there a way to convert LAS point clouds to DWG and/or DXF files using open source software?

Also willing to script with Python - is this possible using libLAS/LASlib for example?


texts: I have the OpenLayers plugin on my computer (user/.../plugins), but it does not show when I open the plugins dialogue box. How do I get it to show in the plugins menu?


texts:
I have a custom map, served by our custom internal map server (XYZ) showing up in OL3. This works perfectly, but now I have to put some vector shapes over the top, these shapes are coming from an SVG. Both the tiles and the shapes which are being overlayed are coming from the same source SVG (map server tiles converted to PNG at various zoom levels).

The issue that I have is that the shapes (now converted to GeoJSON for OL3) are appearing upside down, the Y-Axis is inverted through the centre of the map. I have proven this by manually inverting the shape in gimp (flip through the middle of the image) and it then appears as it should.

Is there a way to get OL3 to do this for me or do I have to code something?  Both the extent and the SVG shapes have the bottom left corner as 0,0 - the extent therefore is:

var extent = [0, 0, 8192, 8192];

The projection I am using for OL3 is below:

var projection = new ol.proj.Projection({
    code: 'EPSG:4326',
    units: 'pixels',
    extent: extent
});

The shape from the SVG (converted to GeoJSON) is:

{
        "type": "Feature",
        "id": 1006,
        "properties": {
            "extId": "svg6",
            "type": "test shape",
            "name": "svg6",
            "category": "test category"
        },
        "geometry": {
            "type": "Polygon",
            "coordinates": [[[4644.3125,6600.9688],[4643.125,7269.5938],[4577.4688,7270.9688],[4578.6562,7990],[5740.5625,7990],[5740.5625,8177.9062],[6298.625,8179.2812],[6298.625,7833.375],[6055.9375,7834.75],[6057.125,6609.125],[4644.3125,6600.9688]]]
        }
    }

Hope that makes sense, let me know if anyone needs any more info...


texts: How to display the geometry type of features in the attribute table?

texts: Does anybody know how to display the geoemtry type of features in the attribute table?


texts: Calculating mean raster value within small polygon zones using ArcGIS Desktop

texts:
I am trying to calculate the mean raster value within zones (vector) using the tool "zonal statistics as table". This works for the larger-area zones, but not for smaller-areas. Making the cell-size of the IDW smaller does work.

How can I determine what cell size to use?

But, I am looping through many files so would prefer to not to do this since it takes much longer.

Are there alternative tools?


texts:
I am using the following code to extract raster values over a SpatialPolygonDataFrame.

result<-extract(Stack, polygons)

>Stack
class       : RasterStack
dimensions  : 634, 862, 546508, 13  (nrow, ncol, ncell, nlayers)
resolution  : 10, 10  (x, y)
extent      : 492660, 501280, 4378960, 4385300  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : L2A_T29SM//21_B02_10m//etc.
min values  :             1.0000000,             1.0000000,             1.0000000,             1.0000000,             1.0000000,             1.0000000,             1.0000000,           119.0000000,            88.0000000,             1.0000000,   -0.9932203,   -0.9959184,    1.0000000
max values  :          1.460000e+04,          1.547100e+04,          1.485900e+04,          1.547600e+04,          1.168700e+04,          1.131400e+04,          1.109100e+04,          1.084500e+04,          1.451000e+04,          1.105000e+04, 9.996238e-01, 9.992918e-01, 4.000000e+00

> polygons
class       : SpatialPolygonsDataFrame
features    : 13
extent      : 493428.6, 499603.9, 4379481, 4383781  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=29 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
variables   : 3
names       : id, ID_1,    LC
min values  :  1,    1,  Agri
max values  :  9,    4, Water

Then I want to create a dataframe out of it but I received the following error:

resultdf<-as.data.frame(result)
Error: arguments imply differing number of rows: 225, 320, 286, 351, 68, 79, 41, 63, 70, 171, 72, 158, 497.

From what I check in the forum it is a problem with the data structure since the dataframe would not have the same length; how can I solve it?


texts: Can't import certain tables from SQL Server to QGIS

texts:
This is a list of 193 zip codes and their geometries that do not import in QGIS.

This is a list of 999 Census Block Groups and their geometries that do import in QGIS.

This is the Add MSSQL Layer tbl in QGIS. It shows what the table displays when I try to import them both.

What would cause one table to be able to import, but not the other?

I am having this problem where I can't import shapefiles from SQL to QGIS. I am using zip codes in Florida as my shapes. I wrote a query in sql to give me certain zip codes and their geometries for certain counties in FL. This was put into a new table in sql consisting of nearly 200 zip codes. When I try to import this into QGIS, it wont import. The MSSQL import table doesnt allow me to select this table. However, if I change the query to include only 20 zip codes, it will allow me to select and import the table. So does QGIS have a limit to the amount of shapes it can import? I have other tables in sql that have more than 1000 shapes and they import, but not this one of 200 zip codes.


texts: Convert coordinates with Proj4Leaflet from EPSG:25832 to leaflet 3857

texts:
I have a wmts server, and try to add this tiles to leaflet with proj4leaflet.
My code is here

//this one is fine
     var main = L.tileLayer(
          'http://****&zoom={z}&x={x}&y={y}', {
            attribution: 'Maps'
          })

//this one uses EPSG:25832
     var url =
          "https://****/gk/gk.nib_utm32_wmts_v2?&gkt=673B1E78E90DBCB2073CC248B949B95CE5156C8E47F6E53B68C0F2A903F3C4C1BC9F3C9E2E1D16E437B7E784B683185C5A8F6E3ADF7BD0828D018F0853CD0DA1";
        var ign2 = new L.TileLayer.WMTS(url, {
          layer: "UTM32_EUREF89",
          style: "default",
          tilematrixset: "default028mm",
          format: "image/png",
          height: 256,
          width: 256,
        });

        var crs = new L.Proj.CRS('EPSG:25832', "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs",new L.Transformation(1, 0, -1, 0));
        map = new L.Map('mapid', {
          crs: crs,
          continuousWorld: true,
        });
        map.addLayer(main);
        map.addLayer(ign2);
        map.setView([50,4], 4);
        L.marker({
          "lat": 59.87891,
          "lng": 8.59411
        }).bindPopup('Test projection').addTo(map)

So this code does not work. I got an error saying
Uncaught TypeError: Cannot read property '4' of undefined on this line map.setView([50,4], 4);

When I comment this line crs:crs here

 map = new L.Map('mapid', {
              //crs: crs,
              continuousWorld: true,
            });

I got this result, so I have to scale down and move that big Norway.

Got this info in tiles xml

<TileMatrix>
<ows:Identifier>0</ows:Identifier>
<ScaleDenominator>7.737142857141884E7</ScaleDenominator>
<TopLeftCorner>-2000000.0 9045984.0</TopLeftCorner>
<TileWidth>256</TileWidth>
<TileHeight>256</TileHeight>
<MatrixWidth>2</MatrixWidth>
<MatrixHeight>2</MatrixHeight>
</TileMatrix>
<TileMatrix>
<ows:Identifier>1</ows:Identifier>
<ScaleDenominator>3.868571428570942E7</ScaleDenominator>
<TopLeftCorner>-2000000.0 9045984.0</TopLeftCorner>
<TileWidth>256</TileWidth>
<TileHeight>256</TileHeight>
<MatrixWidth>3</MatrixWidth>
<MatrixHeight>3</MatrixHeight>
</TileMatrix>
<TileMatrix>

UPD:
Added resolutions from xml ScaleDenominator

var resol = [77371428.57141884,38685714.28570942,19342857.14285471,9671428.571427355,4835714.285713677,2417857.1428568386,1208928.5714284193,604464.2857142097,302232.14285710483,151116.07142855242,75558.03571427621,37779.017857138104,18889.508928569052,9444.754464284526,4722.377232142263,2361.1886160711315,1180.5943080355657,590.2971540177829];
    var crs = new L.Proj.CRS('EPSG:25832', "+proj=utm +zone=32 +ellps=GRS80 +units=m +no_defs", {
      resolutions: resol
    });
    map = new L.Map('mapid', {
      crs: crs,
      continuousWorld: true,
    });

Nothing( 400 (Bad Request) from server). But with resol.reverse() got this, and zoom out when zoom in (that is because of reverse as I understand)


texts:
I would like to use a set of color ramps that are available at CPT-CITY but cannot find any way to import the downloaded files.
Note that I have found the Catalog: cpt-city option under the Create New Color Ramp selection. However, the ramps that I'm interested are not found in this collection.
Is there some way to simply browse to the file that I'm interested in to register it with my QGIS project?


texts: Append/Insert geoprocessing results into an existing feature class (polyline midpoints)

texts:
I plan to run a scheduled server job on a feature class in an Oracle 18c/10.7.1 geodatabase.

The job would:

Truncate the rows in an existing feature class
Generate new rows from a geoprocessing tool (FeatureVerticesToPoints / midpoints)
Insert/append the geoprocessing rows into the feature class

Note: I would like to do this without needing to temporarily store any sort of intermediate results in a temporary feature class/GDB.

I've gotten part-way there:

import arcpy
conn = "Database Connections\\DEV.sde\\"
input_lines = "OWNER.ROAD"
existing_fc = "OWNER.SCHED_ROAD_MIDPOINT"

arcpy.TruncateTable_management(conn + existing_fc)
arcpy.FeatureVerticesToPoints_management(conn + input_lines, conn + existing_fc, "MID")

#arcpy.Append_management(, , "TEST", "", "")

Question:

I'm not sure how to use FeatureVerticesToPoints_management as the input in the append process.

Is there a way to append/insert geoprocessing results into an existing feature class using ArcPy?


texts: Some related table fields not displaying in ArcPad 8

texts:
I have an issue that I have been struggling with and I'll try to explain without getting too confusing.  I have an AXF that is displaying only some of the fields in a related table on the default attributes page for related tables.  So here are the details, I have 6 fields in the related table and a relationship class that relates the table to the relationship class.  Here are the specifics of the relationship;  feature class is the origin, related table is the destination, relationship is composite, messages are forward, 1-M, Primary key is GlobalID (GUID) to AuditID (GUID).  The AXF is checked out from ArcMap.  When I open the AXF in ArcPad 8, and select a feature to view the related table attribute page, the list view page is only displaying 4 of the 6 fields (Inspector - text field, Pass/Fail - coded domain, contractor - coded domain, and comments - text field).  When I switch to single line view to add data, the only fields I see are Inspector, Pass/Fail and contractor.  The comments field is not displayed (even though it was displayed on the list view page).  The other fields in the table (AuditID - GUID primary key, AuditDate - date) are never displayed on either view of the attributes table.  When I first developed the related tables and relationships, everything was working as expected with related table fields displaying in ArcPad.  When I created the production version of the geodatabase and added the related tables and relationships, things are not displaying the same.  Any yes, I created the prod database the same way I create the development one.  If any one has any insight or an answer to what could be happening, I would really appreciate the help.  My client is waiting for this application and I need to get out to the customer ASAP.

Thanks.

Tim


texts: Openlayer3 3.9 vector layer ignores some linestrings

texts:
I updated my small OL3 project to 3.9 but noticed that one of my vectorlayers now only displays very few features:

Other point layers work fine. Also I checked the basics:

Browser calls and transfers geoJSON successfully
my geoJSON is valid and get displayed at geojsonlint.com (using leaflet)
firebug doesn't fires errors, even with ol-debug

As said it worked fine with OL3.5 but I needed to add the GeoJSON format to migrate. Works fine for point layers :/
Here is my code snippet:

    ...
    ,
    new ol.layer.Vector({
        title: 'Links',
        source: new ol.source.Vector({
            url: '/api/links',
            format: new ol.format.GeoJSON({
                //defaultDataProjection :'EPSG:4326',
                projection: 'EPSG:3857'
            })
        }),
        style: [new ol.style.Style({
            stroke: new ol.style.Stroke({
              color: 'yellow',
              width: 5
            })
          })],
          visible:true
    }),
    ...

Has anybody ideas how to isolate the problem? I guess I didn't found an OL bug, did I?


texts: Obtaining coordinates from center of grid squares (research tools)

texts: I have a grid from research tools saved as shapefile. I want to display the center coordinate of each grid square in the attribute table.


texts: 3D rubber sheeting using ground control points

texts:
I have 3d lines that I would like to rubber sheet or rectify in 3D (including z rectification). I used ogr2ogr with -gcp, however all the vertex of a polyline were set to a single elevation, and since you don't input a z for reference I am guessing it does not rubber sheet in 3D.

ogr2ogr feat_rect.shp feat_in.shp -gcp 765825.824 2169178.054 765820.33 2169177.15 413.5 -gcp 765849.205 2169176.512 765846.487 2169174.837 -gcp ...

It seems like GRASS v.rectify could do it, however I'd like a solution where I don't need to setup a GRASS workspace and rather run it in batch mode. ESRI solution could also be acceptable if I can run it in batch mode.


texts: Trying to write standalone Python script for QGIS, can't load SIP

texts:
I'm trying to learn to write standalone Python scrips for QGIS, but can't get past first base.

I'm using QGIS 2.18, installed with the OSgeo4W64 installer.  I have modified my path as follows.

SET OSGEO4W64_ROOT=C:\OSGeo4W64
SET QGIS_PREFIX=%OSGEO4W64_ROOT%\apps\qgis
SET PATH=%PATH%;%QGIS_PREFIX%\bin
SET PYTHONPATH=%QGIS_PREFIX%\python;%PYTHONPATH%

I'm trying to run a script that begins:

import os, os.path, sys
from qgis.core import *
from qgis.gui import *
from PyQt4.QtGui import *
from PyQt4.QtCore import *

I get the following error:

C:\Users\Steve\Documents\MEGAsync\Steve\Lex>lex.py
Traceback (most recent call last):
  File "C:\Users\Steve\Documents\MEGAsync\Steve\Lex\lex.py", line 3, in <module>
    from qgis.core import *
  File "C:\OSGeo4W64\apps\qgis\python\qgis\__init__.py", line 26, in <module>
    from qgis.PyQt import QtCore
  File "C:\OSGeo4W64\apps\qgis\python\qgis\PyQt\QtCore.py", line 26, in <module>
    import sip
ImportError: DLL load failed: %1 is not a valid Win32 application.

What am I doing wrong?


texts: Spatial join attributes in QGIS for a time series

texts:
I have metropolitan areas and I want to calculate average values (GDP, Population, etc) per year. I have another layer with NUTS areas with theses values, from 1970 to 2012.

So far, I know the operation "join attributes by location" to calculate for example the average of GDP per inhabitant or the sum of the population, but I have to repeat the operation for each year.

Is there a way to do it automatically ? Or even better, to get the percentage of the area for each NUTS polygon layer ? Like metropolitan area MA1 is 40% NUTS N1, 55% N2 and 5% N3, so I would be able to calculate the averages and the sums.

I am a beginner with QGIS. I apologize if the question has been already answered, I have not found anything so far about that.


texts: What would be a good way to test data in different projections have identical locations?

texts:
I have 2 sets of files for which I want to compare the coordinates of ~3000 plots

I have a 38 shapefiles of surveyed plots in different regions of Northern Canada. Some of them are projected in Lambert Conformal Conic, others are projected in Albers Equal Area.
I have an Excel database with these same plots. However, in this database, the plots are provided in UTM. They cover 13 different UTM zones

I want to know if the locations of the plots are the same or different between the shapefiles and the Excel file. I'm expecting to find instances of both, and I want to know for which plots the coordinates have changed.

My solution to try to figure this out was to use R to measure the distance between each plot's location in the shapefiles and the Excel file. If the result is 0 , the location is correct. If the result is >0m, I need to identify which location is correct, and change the one that is incorrect. Realistically, I could probably count them as equivalent if they're less than 10-15m off.

I had individually converted the UTMs from each zone to Albers Equal Area, and then all of the shapefiles to Albers Equal Area. However, someone mentioned to be that Albers Equal Area is good for measuring area, but not distance.

I thought I would re-run my code, but convert everything to Lambert Conformal Conic instead. However, read a source saying: "Conformal projections preserve the correct angles between directions within small areas, though distorting distances."

I'm stumped. What WOULD be a good projection for measuring distances?

Or alternatively, is there a better way to approach the problem of testing if  data in different projections have identical locations?

(I would prefer to use R. Usually I have access to ArcGIS, but I don't right now due to the pandemic. I have downloaded QGIS but I'm not very comfortable with it)


texts:
I have created a kml file with around 300 place markers.  When I try to display the kml by typing the address in the maps search bar, the map seems to behave erratically:
- If I click on the marker sometimes it shows a blank callout.  If I click on the marker in the side bar it shows the bubble correctly. If I zoom in and click on the marker in the side bar, the callout appears in a random location.

Is this because of too many markers?

Here is my file:
https://sites.google.com/site/nilanjanmaps/home/testmap/softwaretesters2012all.kml


texts:
I know nothing about GIS, but I have what I think is a GIS problem.  From spending a few hours googling my problem, I think what I need to do is perform a spatial query, but I'm not sure how to go about it.

I have some crime data (a list of specific incidents, the time date that they occurred and the address where they occurred) that I'd like to look at in relation to one particular neighborhood. I have an outline of the neighborhood on a map, which is specific enough that I could translate it into a polygon that describes the outline of the neighborhood.  I'd like to filter the crime data I have to see which specific incidents happened within the bounds of this neighborhood.  I think this is a spatial query, but I'm not sure.

I'm looking for some guidance on the problem in general.  It appears that there are some free tools online for these kinds of things (Google Fusion Tables?) and I have access to ArcGIS Desktop and MS SQL Server.  I'm familiar with SQL Server, but not ArcGIS Desktop and it doesn't look like something I could pick up and learn in a day or two, but I'm not sure.


texts: sf::st_transform() gives "crs not found" error

texts:
I am using sf package version 0.9-6. I define a new object in ESRI:4326 and then try to reproject it to an Albers projection. I get an error, despite explicitly setting the crs when defining the object.
Here is MWE.
library(sf)
# Linking to GEOS 3.8.1, GDAL 3.1.4, PROJ 6.3.1
disp_win_wgs84 <- sf::st_sf(st_sfc(st_point(c(-97.2, 32.55)),
                         st_point(c(-85.55, 49.2)),
                         crs = 4326))
target_crs <- 102003
disp_win_trans <- sf::st_transform(disp_win_wgs84, crs = target_crs)

I get the following error:
Error in CPL_transform(x, crs, aoi, pipeline, reverse) :
  crs not found: is it missing?
In addition: Warning message:
In CPL_crs_from_input(x) :
  GDAL Error 1: PROJ: proj_create_from_database: crs not found

Edit: updated to define the target_crs as a number. Apparently, this is what was causing the error (see the answer). The solution is to use target_crs <- "ESRI:102003", this does not raise an error.


texts: Any source code for converting an infix filter to a WFS filter?

texts:
A WFS GetFeature query looks something like this:

http://hname/arcgis/services/sname/MapServer/WFSServer?SERVICE=WFS&REQUEST=GetFeature&VERSION=1.0.0&MAXFEATURES=1001&FILTER=%3CFilter%3E%3CPropertyIsEqualTo%3E%3CPropertyName%3EADDRKEY%3C/PropertyName%3E%3CLiteral%3E1001%3C/Literal%3E%3C/PropertyIsEqualTo%3E%3C/Filter%3E&TYPENAME=Addresses

But users would want to enter something like "ADDRKEY = 1001" so it becomes necessary to convert the infix expression to a AST and then into a WFS filter.

Are there any reusable solutions for converting an CQL filter into a OGC filter?

"A=b" is converted into
"<PropertyIsEqualTo><PropertyName>A</PropertyName><Literal>b</Literal></PropertyIsEqualTo>"


texts:
I am new to GIS and I want to do Suitability Analysis for Landfill Site Selection in QGIS. I had a Land-Use map in jpeg format, from which I started by making shapefiles of each kind of input layer such as commercial areas, settlements, agriculture, roads, water bodies etc. Now, I am having problems in rasterizing these shapefiles. I have not yet created any fields in the attribute table as I do not know exactly what is the role of fields in rasterizing.
I tried to rasterize by using Raster menu> Conversion> Rasterize. but it produced just a black rectangular layer with nothing on it. Here are the Screenshots of one of the shapefiles and the raster that got created by following above method-


texts:
I have the following problem in QGIS:

I try to create voronoi polygons based on a shapefile (points; standard WGS84). The newly created voronoi polygons layer seems to be ok for the WGS84 CRS. But when I swith into Pseudo Mercator CRS (3857) then the voronoi polygons are (of course) not representing the correct voronoi polygon layer and the polygons are not representing areas according to the neares points of the souce layer for that CRS.

Is there a way, how to create new voronoi polygons layer based on the same source shapefile (standard WGS84), but so that the final voronoi layer would work correctly in the Pseudo Mercator projection?


texts: How to convert jpgs to RST for IDRISI Selva Image Processing

texts: I would like to create a mosaic data set in IDRSI Selva however, my images are jpg.  Is there a way to convert these within IDRISI or otherwise?


texts:
I have two shapefiles. One is a point layer of trees, the other is a polygon of cadastral boundaries. I want to spatially join these two layers so that I get the information, on which cadastral boundary the tree lies. For every point of trees, I want to add the information of the boundaries to the trees in the query.

I have seen something similar here, but I'm not a big SQL expert :-)

Spatial query, embed in layer or project

Perhaps somebody can help me out.


texts:
I want to open a csv raster data file with GDAL in C#

Dataset csvFile = Gdal.Open(fileLocation, Access.GA_ReadOnly);

Produce an error saying that

At line 4, X is 412267, where as 412555 was expected

I tried to change the X manually but instead it is saying

At line 4, X is 412555, where as 412555 was expected

What confused me more is when I tried to open this csv file with C++ or python it works fine.

Here the python code I used

dataset = gdal.Open( filename, GA_ReadOnly )

Anyone here facing the similar problem with me?


texts: Map box GL mature / appropriate for own OSM-Tileserver

texts:
I am quite new to this topic and currently trying to decide what library to use for developing OSM-apps with my own tile-server for Android, iOS and Websites (currently using Leaflet for websites) - Windows Phone might follow, but that's on another sheet.

I was impressed by the features of Mapbox GL - but I am not sure if it is ready for productive use.

1) Is this stable and performant or should I use another library (e.g. MapBox)?

2) I am also not totally clear if I can use MapBox or MapBox GL with my (automatically updating) own OSM-Tileserver without purchasing a plan from MapBox. My current understanding what could work is :

Download OSM XML
Convert it with osm2pgsql (convert data from osm to osm tiles/marble)
Convert it with TileMill API (How to automate export in Tilemill?)

Is this the way to go?


texts: How can I get nearest road to a lat lon in MSSQL?

texts:
I've updated osm data (ways) in MSSQL and got geographic objects. Now I want to get the closet road to any lat lon a user will query.

Which SQL spatial statement is used for that ?


texts: Count unique points in polygon Python Error

texts:
I need to count different values in a point shapefile field inside several thousands of polygons. I know the tool "Count unique points in polygons" and it works fine when I call it from the toolbox. But as I said, I have to use it in a lot of polygons, so I wrote this Python script:

import os, processing

# Set local variables
path = r'my\base\path\Dissolves'
continents = r'\somepath\points.shp'
outpath = r'path\to\outfile'

# Creamos la lista de shapes
shpList = [os.path.join(path, i) for i in os.listdir(path) if i.endswith('.shp') and not 'point' in i]
print('creada lista de', len(shpList), 'shapes')

for shp in shpList:

    outShp = os.path.join(outpath, os.path.split(shp)[1])
    print('Haciendo especie:', os.path.split(shp)[1])

    try:

        # Ejecutamos el join espacial dentro de un tr por si hubiera algun error poder continuar al siguiente shape
        processing.runalg('qgis:countuniquepointsinpolygon', shp, continents, 'issed',  'NUMPOINTS', outShp)

    except Exception as e:

        print('Error', e, os.path.split(shp)[1])
        continue

When I applied this code I don't get any error message, and I get the new polygons and everything seems to be ok. But I don't know why the higher value in the count unique values field ("NUMPOINTS") it never gets higher than 2, which is incorrect.

I tried the same polygons with the same points file, doing them one by one and doing them in the script, and get different values (I mean, all the values larger than 2 appears to be 2 in the script).

I edit this post to write the link to some shapes (just 5 and a part of the points in my Dropbox: https://www.dropbox.com/sh/f5982yd8dovsng3/AAAfwVmbtGc4uSOhi96GLHP4a?dl=0).


texts: Starting ArcMap editing session without displaying Editor toolbar?

texts:
I'm developing a toolbar with some geoprocessing routines and digitizing functionality. On this toolbar I have implemented a "start edit" button.

I want to start an edit session by clicking that button on my toolbar, but the original ArcMap Editor toolbar shows up too, even so I deselected it form the customize menu.

Is it possible to prevent the Editor toolbar from opening, so users are forced to use my custom toolbar?

You can simulate my workflow by adding the "Start Editing" button to any toolbar in ArcMap's Customize mode. Note that the Editor toolbar always switches on when starting an edit session.


texts: Reduce number of null values in MODIS surface temperature data

texts:
I am trying to reduce the number of null (no data) values in output data for an area (Tokyo prefecture) in August 2019 when using MODIS 1km spatial resolution daily surface temperature satellite data (MODIS/006/MOD11A1, LST_Day_1km).

I am new to Google Earth Engine and coding but through a combination of watching tutorials, trial and error and cut and pasting from code structures found online (especially from a question asked on StackExchange by GGamer and answered by Kuik) I have managed to produce the following code;

https://code.earthengine.google.com/7e467c91855832cce3b65c5dfce4dd1e

However my output data (in degrees celsius) for each point on each day in August in the Tokyo area always has a large number of nulls, for example;

FeatureCollection (2100 elements, 32 columns)
type: FeatureCollection
columns: Object (32 properties)
features: List (2100 elements)
1: Feature 1 (Point, 31 properties)
type: Feature
id: 1
geometry: Point (139.01, 35.89)
properties: Object (31 properties)
08_01_2019-08-01: null
08_02_2019-08-02: 24.670000000000016
08_03_2019-08-03: null
08_04_2019-08-04: 24.55000000000001
08_05_2019-08-05: 24.250000000000057
08_06_2019-08-06: 26.130000000000052
08_07_2019-08-07: 24.910000000000025
08_08_2019-08-08: null
08_09_2019-08-09: 26.57000000000005
08_10_2019-08-10: null
08_11_2019-08-11: null
08_12_2019-08-12: null
08_13_2019-08-13: null
08_14_2019-08-14: null
08_15_2019-08-15: null
08_16_2019-08-16: null
08_17_2019-08-17: null
08_18_2019-08-18: null
08_19_2019-08-19: null
08_20_2019-08-20: null
08_21_2019-08-21: null
08_22_2019-08-22: null
08_23_2019-08-23: null
08_24_2019-08-24: 18.010000000000048
08_25_2019-08-25: null
08_26_2019-08-26: null
08_27_2019-08-27: null
08_28_2019-08-28: null
08_29_2019-08-29: null
08_30_2019-08-30: null
ID: 1

Is there something I can add to/remove from my code to reduce the amount of nulls? Or is this unavoidable? Or perhaps combining with a different dataset is another option to fill the gaps, for example MODIS/006/MYD11A1, LST_Day_1km?


texts: GeoServer - Multiple symbology on single layer

texts:
I'm trying to show multiple symbologies on a single GeoServer layer.
I'm also using the CSS extension, but an SLD example would be fine.

My GeoServer version is 2.8.2

For example:
A yellow rectangle with a red circle within it(on a point layer).

Is there a way for me to achieve this?


texts: GeoDjango web mapping application can only view single point at one time?

texts:
I am very new to Geodjango and have been running into some issues with a web mapping application. Currently, I am trying to make it so I can view all the schools on a map at once, as well as view individual schools on a map. I have been able to make it so I can view individual schools, but I am still unable to view all schools on a map at once.

This is how it looks when I am able to view a single school:

Here is what happens when I try to view all schools at once:

Am I going to the wrong url to do so or is there some other issue in the code that I am not seeing?

Here are my scripts:

schools_project/urls.py:

from django.conf.urls import url, include
from django.contrib import admin
from django.urls import include, path
...
urlpatterns = [
    url(r'^schools/', include('schools.urls')),
    path('admin/', admin.site.urls),
]

schools/models.py:

from django.db import models
from django.contrib.gis.db import models as geomodels
...
class School(models.Model):
    name = models.CharField(max_length=300, blank=False)
    geometry = geomodels.PointField()
...
class Meta:
    order = ('name',)

schools/apps.py:

from django.apps import AppConfig
...
class SchoolsConfig(AppConfig):
    name = 'schools'

schools/admin.py:

from django.contrib import admin
from .models import School
from leaflet.admin import LeafletGeoAdmin
...
class SchoolAdmin(LeafletGeoAdmin):
    list_display = ('name', 'geometry')
...
admin.site.register(School, SchoolAdmin)

schools/urls.py:

from django.conf.urls import url
from . import views
...
app_name = 'schools'
...
urlpatterns = [
    url(r^school/(?P<pk>[0-9]+)$',
        views.SchoolsDetailView.as_view(), name='school-detail'),
]

schools/views.py:

from django.views.generic import DetailView
from .models import School
...
class SchoolsDetailView(DetailView):
    """
        School detail view.
    """
template_name = r"C:\Users\Drew\Desktop\djangoproject\schools_project\templates\school-detail"
    model = School

school-detail.html (template):

<!DOCTYPE html>

<Title>Public Schools in Aurora, Colorado</title>
<meta name="viewport" content="width=device-width, initial-scale=1">
<link rel="stylesheet" href="https://www.w3schools.com/w3css/4/w3.css">
{% load leaflet_tags %}
{% leaflet_css %}
{% leaflet_js %}

Public Schools in Denver and Aurora, Colorado

{{ school.name }}

    {% leaflet_map "main" callback="map_init" %}

Does anybody have any ideas or tips on how I can view all points on a map at the same time?


texts: Using R to create a choropleth generated from lat/lon lists - with data generated internally

texts:
a friend referred me here as opposed to the regular Stack Overflow site, as this community maybe better able to help!

My issue is largely with the syntax and framework of ggmap,ggplot, and the sp package.

I have already summarised and prepped my data. This link shows (some of) the data I'm trying to map. I used the following code block to create Spatial polygon objects instead of lists of lat/long values, as shown in the data in the above link:

createShape = function(sub){
  #This funciton takes the list of lat/lng values and returns a SHAPE which should be plottable on ggmap/ggplot
  tempData = as.data.frame(do.call(rbind, as.list(VICshapes[which(VICshapes$Suburb==sub),] %>% select(coords))[[1]][[1]]))
  names(tempData) = c('lat', 'lng')
  p = Polygon(tempData)
  ps = Polygons(list(p),1)
  sps = SpatialPolygons(list(ps))
  return(sps)
}
VICshapes[['shape']] = lapply(suburbNames, function(x) createShape(x))

So - I've got the data I want to plot. I've got shapes. I've got the packages installed... but how do I actually make things look pretty and presentable? ANY of the help on packages seems to assume I have .shp files, along with other external data which defines how I plot it all. Apparently I need to use ggmap::fortify() or broom::tidy on the shape files to get them to plot. Then I just get warnings from ggmap and ggplot about the length of aes arguments being incorrect.

I'm genuinely not sure how to proceed with this, but I know I have the necessary data to show what I want to show - I just don't know how to move forward.

Using R's inbuilt plot function on the output of the createShape() function above shows me a shape as I want. Great - but what if I want ALL the shapes, then coloured according to another variable - which is linked by a common key?

At this stage I'm so frustrated with it all - I'll literally pay someone's rent for a month if you can help.


texts:
I was wondering if there is a LayerSwitcher control in OpenLayers 3?

I read that ol3 is a complete rewrite. But I can't seem to find the equivalent of LayerSwitcher from http://openlayers.org/en/latest/apidoc/ol.control.html. I don't see how one can turn on and off layers without this control. Is LayerSwitcher renamed to something else, missing or is the documentation incomplete?

If OpenLayers3 does not have LayerSwitcher, does anyone know of an example that implements a custom control like the old LayerSwitcher?


texts:
In this script, I take an image and mask off everything below certain elevation.  The result looks really interesting as one can still see the names of the towns and roads where the image is masked off.

var geometry =
    /* color: #98ff00 */
    /* displayProperties: [
      {
        "type": "rectangle"
      }
    ] */
    ee.Geometry.Polygon(
        [[[-83.87143189191625, 28.313819801664945],
          [-83.87143189191625, 24.92605449144572],
          [-78.66391236066625, 24.92605449144572],
          [-78.66391236066625, 28.313819801664945]]], null, false);

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

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

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

  return image.updateMask(mask).divide(10000);
}

// Map the function over one year of data and take the median.
// Load Sentinel-2 TOA reflectance data.
var dataset = ee.ImageCollection('COPERNICUS/S2')
                  .filterDate('2018-01-01', '2018-06-30')
                  // Pre-filter to get less cloudy granules.
                  .filter(ee.Filter.lt('CLOUDY_PIXEL_PERCENTAGE', 20))
                  .map(maskS2clouds);

var rgbVis = {
  min: 0.0,
  max: 0.3,
  bands: ['B4', 'B3', 'B2'],
};

var coper2 = dataset.median()

Map.centerObject(geometry);
//Map.addLayer(coper2.clip(geometry), rgbVis, 'RGB');

var elevation = ee.Image('USGS/SRTMGL1_003')

var above0 = coper2.updateMask(elevation.gt(0))
var above10 = coper2.updateMask(elevation.gt(10))
var above100 = coper2.updateMask(elevation.gt(100))

//Map.addLayer(elevation.neq(0))
Map.addLayer(geometry, {color:'blue'}, 'inundated')
Map.addLayer(above10.clip(geometry), rgbVis, 'above10')

var img = above10.visualize({
  bands:['B4', 'B3', 'B2'], min: 0, max: 0.3
    })

Export.image.toDrive({
  image: img,
  description: 'florida_above_10',
  scale: 1000,
  maxPixels: 3e9,
  region: geometry})

Here's what is seen in GEE in the browser:

with the geometry added as a map layer, and

with no geometry added as a map layer.

Is it possible to export an image file that will be able to show what's seen in either of these when opened (perhaps with the map in the background)?


texts: Convex hull with n points using computational geometry libraries

texts:
I've been searching for a function to make the Convex Hull of a Spatial Object with N points, to make an spatial approximation with n-corners, as shown in the article "Comparison of Approximations of Complex Objects Used for Approximation-based Query Processing in Spatial Database Systems", by Thomas Brinkhoff, Hans-Peter Kriegel and Ralf Schneider.
I haven't found any function to make this, I searched for any GEOS function but didn't find. Does anybody know what can I use for it?