texts: Defining spatial reference for WMS layer in ArcGIS API for JavaScript?

texts:
I'm just starting with the ArcGIS Javascript API, and I'm having some issues adding layers from my geoserver to my basic map application.

I started with the WMS resource info sample from the ESRI site (https://developers.arcgis.com/javascript/jssamples/layers_wmsresourceinfo.html). I took out the two layers they display in that sample, and I replaced it with a WMS layer from my geoserver. I'm getting this error:

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://mygeoserver.com/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" >   <ServiceException code="InvalidSRS">
  Error occurred decoding the espg code EPSG:102100
No code &quot;EPSG:102100&quot; from authority &quot;European Petroleum Survey Group&quot; found for object of type &quot;IdentifiedObject&quot;.
</ServiceException></ServiceExceptionReport>

The problem seems to stem from the difference in spatial references. The base map uses 102100, and my layer uses 4326. It doesn't look like changing the spatial reference of my layer is an option, so is it possible for me to change the spatial reference of the base map?

I added 102100 as a custom CRS to my geoserver. I then changed the declared SRS for my layer to EPSG:3857, and told it to "reproject native to declared" (I'm also new to using geoserver, so please let me know if it sounds like I'm doing something wrong there). I'm getting a different error message now:

<ServiceExceptionReport xmlns="http://www.opengis.net/ogc" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" version="1.3.0" xsi:schemaLocation="http://www.opengis.net/ogc http://localhost:8040/geoserver/schemas/wms/1.3.0/exceptions_1_3_0.xsd">
<ServiceException code="InvalidCRS">
Error occurred decoding the espg code urn:x-ogc:def:crs:EPSG:102100 Error in "PROJECTION": No transform for classification "Mercator_Auxiliary_Sphere". Error in "PROJECTION": No transform for classification "Mercator_Auxiliary_Sphere". No transform for classification "Mercator_Auxiliary_Sphere".
</ServiceException>
</ServiceExceptionReport>

And for full information, here's the js code of my basic ArcGIS map:

var map;
require([
  'esri/map', 'esri/layers/WMSLayer', 'esri/layers/WMSLayerInfo', 'esri/geometry/Extent',
  'dojo/_base/array', 'dojo/dom', 'dojo/dom-construct', 'dojo/parser',
  'dijit/layout/BorderContainer', 'dijit/layout/ContentPane', 'dojo/domReady!'
], function(Map, WMSLayer, WMSLayerInfo, Extent, array, dom, domConst, parser) {
  parser.parse();

  map = new Map('map', {
    basemap: 'hybrid',
    center: [-96, 37],
    zoom: 4
  });

  var layer1 = new WMSLayerInfo({
    name: 'states',
    title: 'states'
  });

  var resourceInfo = {
    extent: new Extent(-120.77027961360125 ,35.79436009671527, -120.74876891832204, 35.81224289313929, {
      wkid: 3857
    }),
    layerInfos: [layer1]
  };
  var wmsLayer = new WMSLayer('http://localhost:8040/geoserver/topp/wms', {
    resourceInfo: resourceInfo,
    visibleLayers: ['states']
  });
  map.addLayers([wmsLayer]);

  var details = dom.byId('details');
  domConst.place('Layers:', details);
  var ul = domConst.create('ul', null, details);
  array.forEach(wmsLayer.layerInfos, function(layerInfo) {
    domConst.create('li', { innerHTML: layerInfo.title }, ul);
  });
  domConst.place('WMS Version:' + wmsLayer.version + '
', details);
});

Any idea what might be causing the error?


texts: Calculating area of polygon in WGS84 Shapefile using QGIS

texts:
I am using QGIS 2.12.3 Lyon and I want to calculate the area of my polygon in m2. For that I used the tool provided by QGIS, but the result is not as it should be. I got from QGIS 50 ha, while I got from the TopMaps Viewer 20 ha – I know the result from TopMaps Viewer is right, because I have another area with almost the same size, where I already know the area size.

Then I tried to calculate the area by the field calculator and used the formula $area as explained in other questions in the GIS StackExchange (How to calculate polygon areas and perimeters using QGIS? or How to calculate polygon areas in QGIS?).
When I try this, I can’t click “ok” and I don’t know why.

Probably I should change the KBS of my shapefile, but I really have a big confusion with that. I am using Google maps as background (EPSG: 3857), my shapefile which I created uses WGS 84, EPSG 4326 and the shapefiles I downloaded from public places uses EPSG 31467 Gauss-Kruger zone 3. Is that a problem? I tried to get them on the same KBS, but then they are at the wrong place. Can that be the problem why the calculation of my area size is wrong?


texts: Getting Land Use data from National Land Cover dataset

texts:
I'm downloading data from the National Land Cover database using the get_nlcd tool from the FedData package. Here is an example for documentation for this package

library(FedData)

vepPolygon <- polygon_from_extent(raster::extent(672800,740000,4102000,4170000),
                              proj4string='+proj=utm +datum=NAD83 +zone=12')

# Get the NLCD (USA ONLY)
# Returns a raster
Veg <- get_nlcd(template=vepPolygon, label='VEPIIN',force.redo=TRUE)

#Get Land Use Stats

LandUseStats<-as.data.frame(Veg@data@attributes)
LandUseStats<-LandUseStats[!LandUseStats$Count == 0,]

LandUseStats

  ID OID Value      Count Red Green Blue    NLCD.2011.Land.Cover.Class Opacity
1   0   0     0 7854240512   0     0    0                  Unclassified       1
12 11  11    11  469012527   0     0    0                    Open Water       1
13 12  12    12    1599206   0     0    0            Perennial Snow/Ice       1
22 21  21    21  292251633   0     0    0         Developed, Open Space       1
23 22  22    22  131633826   0     0    0      Developed, Low Intensity       1
24 23  23    23   59456652   0     0    0   Developed, Medium Intensity       1
25 24  24    24   21426522   0     0    0     Developed, High Intensity       1
32 31  31    31  110507264   0     0    0                   Barren Land       1
42 41  41    41  973617734   0     0    0              Deciduous Forest       1
43 42  42    42 1037912310   0     0    0              Evergreen Forest       1
44 43  43    43  179845520   0     0    0                  Mixed Forest       1
53 52  52    52 1940362409   0     0    0                   Shrub/Scrub       1
72 71  71    71 1306961628   0     0    0                   Herbaceuous       1
82 81  81    81  597234572   0     0    0                   Hay/Pasture       1
83 82  82    82 1392218141   0     0    0              Cultivated Crops       1
91 90  90    90  347144473   0     0    0                Woody Wetlands       1
96 95  95    95  116679631   0     0    0 Emergent Herbaceuous Wetlands       1

#Compare to plot
plot(Veg)

If I understand the "Count" to be the number of pixels represented in the associated plot, then when comparing the plot to the Stats, the two don't seem to match. Nor do I understand why close to half the data is "unclassified".

Is there a better way to extract the land use data used in the plot?


texts: Is there a publicly available street centerline file that can be used to geocode Westchester addresses?

texts: The publicly available file at NYSGIS Clearinghouse does not have the required fields to geocode Westchester addresses. The file available for cooperating members does have them, but my institution is not a cooperative member.


texts: How could you measure the minimum and maximum distance to the outer edge of a polygon in QGIS

texts:
I am looking for a automatic method to measure the minimum and maximum distance from a point (represented by the star) to the outer edge of a polygon (black dotted line).

I know this could be done with the measure tool but this can have accuracy issues and is a manual process.

Are there any plugins that do this of the shelf or could anyone recommend a workflow model or script?


texts: QGIS procedure for Multidirectional, oblique-weighted, shaded-relief?

texts:
Is QGIS capable of the multidirectional, oblique-weighted,shaded-relief method described at  http://pubs.usgs.gov/of/1992/of92-422/of92-422.pdf ?  If so, what are the steps necessary in QGIS to achieve this style of hillshading?

I am able to utilize the hillshade tool in QGIS 2.6 to obtain a "single-directional" hillshade image from 1/3 arc NED data, but I'm interested in the method described in the link.


texts:
I've defined the following map:

Proj4js.defs["EPSG:31258"] = "+proj=tmerc +lat_0=0 +lon_0=13.33333333333333 +k=1 +x_0=450000 +y_0=-5000000 +ellps=bessel +towgs84=577.326,90.129,463.919,5.137,1.474,5.297,2.4232 +units=m +no_defs";
var googleProj = new OpenLayers.Projection("EPSG:900913");
var salzburgProj = new OpenLayers.Projection("EPSG:31258");
var bounds = new OpenLayers.Bounds(426249,294999,426251,295001);
var maxZoom = 19;
var minZoom = 15;
var res = [4.77731426696777,2.38865713348389,1.194328566741945,0.5971642833709725];

mapCOD = new OpenLayers.Map({
    div: "mapCODs",
    allOverlays: true,
    maxExtent: bounds,
    projection: googleProj,
    allowSelection: true,
    controls: [new OpenLayers.Control.PanZoomBar(), new OpenLayers.Control.Navigation()]
    });

var layer = new OpenLayers.Layer.OSM(
            "layer1",
            "tiles/input/${z}/${x}/${y}.png",
            {zoomLevels: maxZoom, zoomOffset: minZoom, resolutions: res}
        );
mapCOD.addLayer(layer);

I have created map tiles by myself using mapnik.
The map is visualized properly, but I have a shift problem when I try to add a vector layer as follows:

// feat = ...well formed GeoJson feature
var geojson_format = new OpenLayers.Format.GeoJSON({internalProjection: googleProj, externalProjection: salzburgProj});
var vector_layer = new OpenLayers.Layer.Vector(1);
vector_layer.addFeatures(geojson_format.parseFeature(feat));
map.addLayer(vector_layer);

You can see the result in the image above: the vector layer (purple geometries) should be up-right shifted, in order to match the tiles below them.
I do not understand why I obtain this behaviour: what is missing in my code?


texts: Keep tracking alive for trackpaths when phone is in sleep mode

texts:
I've found that tracking only works if QField is active and the screen is turned on constantly. This drains battery power. As soon as the power button is pressed to conserve the battery the tracking stops and connects a straight line to the next location when QField is activated again. I am using a bluetooth GPS and have tried Lefebure Design's NTRIP Client as the mock location provider.. despite this you still have to turn on location services for QField to pick up a location. MapitGIS is able to use the mock location provider and it's own GNSS settings to the bluetooth device while recording tracks whether the phone is in sleep mode or not.. this is how most GPS tracking devices work normally (eg Garmin eTrex devices). What is promising about QField's tracking feature is that you can still take points while it keeps tracking.
a) Is there a way to keep the phone's location services turned off but to pick up the location from the mock location provider (eg from Lefebure Design's NTRIP Client)? With location services turned on and airplane mode turned on it does seem to be using the bluetooth device in some way but it theoretically should be able to simply use the mock location provider and leave location services turned off completely to ensure it using only the external bluetooth GNSS device.
b) Is there a way to keep trackpaths tracking when activating sleep mode to conserve power on the device? We typically are out for 8-10 hours surveying and need the trackpaths capturing data every 3 seconds.


texts: QGIS Atlas zoom to intersection between two different layers

texts:
I have two polygon layers, one is about different land owners boundaries and is used for atlas generation, and the other is about forests inside each property. In atlas generation sometimes i want to zoom to the intersection area between the boundary layer (used for atlas iteration) with the forest layer. That is because some times a huge property can have a small part of forest so i want atlas to zoom in to that part of the property. I was thinking to type a function in the scale bar, maybe $atlasgeometry, geometry, intersection, etc. Can anyone help with this issue?

I am using QGIS 2.8.1 Wien.


texts: Clipping raster with another raster, then mosaic both of them?

texts:
I have two raster tiff files one is a raster that I created from 2 Foot contour lines, the other raster file is a small area that was surveyed to half a foot accuracy.

Is there a way to  essentially clip the area out of the 2ft raster where my survey raster resides and then mosaic or merge the survey raster with the 2ft raster?


texts: Error when creating insert cursor in ArcPy?

texts:
I've written several scripts that copy data from MS Access db into a new, empty ESRI Geodatabase.  I've pretty much copied the code and changed it for each table -- the first two scripts have worked as expected, but the third is throwing an error on the line used to create the Insert Cursor.

All of the variables are storing correctly at the top of the script, but the Insert Cursor is never created.

Error Message:

  Runtime error  Traceback (most recent call last):

  File "", line 49, in    File "c:\program files
  (x86)\arcgis\desktop10.1\arcpy\arcpy__init__.py", line 1005, in InsertCursor
      return gp.insertCursor(dataset, spatial_reference)

  File "c:\program files
  (x86)\arcgis\desktop10.1\arcpy\arcpy\geoprocessing_base.py", line
  383, in insertCursor
      self._gp.InsertCursor(*gp_fixargs(args, True)))

  RuntimeError: ERROR 999999: Error executing function.

Code:

#Import modules, environment settings
import arcpy
arcpy.env.overwriteOutput = True

# Database connection variables
geodatabase = 'P:/GIS_Development/ArcGIS Database Design/Landowner.gdb/Parcels'
accessDatabase = 'Y:/Database/LANDOWNER/NewLandowner.mdb/TblParcels'

#Create cursor for ACCESS Database and New Landowner GDB
accessCursor = arcpy.SearchCursor(accessDatabase)

# Create a variable that stores the value for column in a given row in the access table

for accessRow in accessCursor:

    parcelID   = accessRow.getValue('Parcel_ID')
    flag       = accessRow.getValue('Flag')
    agreeID    = accessRow.getValue('Agreement_ID')
    agreeStat  = accessRow.getValue('AgreementStatus')
    ownerType  = accessRow.getValue('Owner_Type')
    add1       = accessRow.getValue('Parcel_Address1')
    city       = accessRow.getValue('Parcel_City')
    county     = accessRow.getValue('Parcel_County')
    state      = accessRow.getValue('Parcel_State')
    zipCode    = accessRow.getValue('Parcel_Zip')
    acreage    = accessRow.getValue('Acreage')
    landClass  = accessRow.getValue('Land_Class')
    mapGrid    = accessRow.getValue('Map_Grid')
    parcelNum  = accessRow.getValue('Parcel_Num')
    parcelCom  = accessRow.getValue('Parcel_Comments')
    enterRes   = accessRow.getValue('Entry_Restrict')
    enterCom   = accessRow.getValue('Entry_Comments')
    contacted  = accessRow.getValue('Contact_Made')
    appStatus  = accessRow.getValue('ApprovalStatus')
    owner      = accessRow.getValue('Owner')

    # "Paste" the stored value from the access table into the new geodatabase
    gdbCursor = arcpy.InsertCursor(geodatabase)
    row = gdbCursor.newRow()

    if state:
        state = state.upper()

    #Fill out fields with variables that are storing MS Access data
    row.Parcel_ID            = parcelID
    row.Flag                 = flag
    row.Agreement_ID         = agreeID
    row.Agreement_Status     = agreeStat
    row.Owner_Type           = ownerType
    row.Address_Line_1       = add1
    row.City                 = city
    row.County               = county
    row.State                = state
    row.Zip_Code             = zipCode
    row.Acreage              = acreage
    row.Land_Classification  = landClass
    row.Map_Grid             = mapGrid
    row.Parcel_Number        = parcelNum
    row.Parcel_Comments      = parcelCom
    row.Entry_Restrictions   = enterRes
    row.Entry_Comments       = enterCom
    row.Contact_Made         = contacted
    row.ApprovalStatus       = appStatus

    # Insert new row into Geodatabase
    gdbCursor.insertRow(row)


texts: How do I convert a fire frequency point shapefile into a fire frequency raster?

texts:
I have a point shapefile that I want to convert to a raster image using ArcMap 10.2. The shapefile has a number of points on top of each other.
Each point is a fire. I want to make a fire frequency raster.

What functions within the conversion toolbox do I use?


texts: How to hide all fields at once in QGIS 1.7

texts:
Is there a way to hide all the fields at once in QGIS 1.7?

I see where I can hide one field at a time but I would like a button which does it all at once.


texts: Ozone concentration when providing coordinates

texts:
Does anybody know where can I find data about ozone? I used TOMS, but it seems that data are not available anymore. What I used to do to retrieve data was specifying lat/lon coordinates and then obtaining the value in Dobson units. I am looking to a service like this one:

http://ozoneaq.gsfc.nasa.gov/ozone_overhead_all_v8.md

But, of course, working!

Thanks


texts: Generating new land cover map using raster and vector data in ArcGIS for Desktop?

texts:
I currently have a fairly large raster (1000, 1000). I have a vector layer of reservoirs. I need to create a 100 m buffer around the reservoir, assign it a land use code and merge it with the larger raster.

I have created the buffer, converted the buffer to raster and then classified it using my land use codes. I then used the mosaic tool to merge the two raster layers. I end up with something like this. The raster layer for my buffer, is a large square and doesn't cover the area of the buffer.

I am wondering if this is because my original raster is so large and the buffer is so small?

What is the best way to go about, creating the buffer, assigning the land use code and merging it with the larger raster?

I am using ArcGIS for Desktop and have access to Spatial Analyst.


texts: Can't access GRASS GIS modules from within Python 2.7

texts:
Using the tutorials posted here, here and here, I am trying to setup a script that will make use of several GRASS GIS raster functions externally. For now I am just trying to get the tutorial script to run but I receive the following error:

CalledModuleError: Module run None ['g.gisenv', '-n'] ended with error Process ended with non-zero return code -1073741515. See errors in the (error) output.

When searching for this error it seems that others have issues trying to import grass.script. However, I am able to import this without issue. It is only afterwards when trying to run a module that I get this error. This happens when I try to use a 'gscript.run_command('g.gisenv') (see code block below).

I run:
- Windows 10 x64
- Python 2.7.13 x86
- GRASS GIS 7.2.0 x86
- x86 version of GDAL

I removed all other GRASS GIS installations (one was installed with QGIS) and reinstalled it with the required Microsoft Visual C++ Redistributable Packages.

import os
import sys
import subprocess

grass7bin = r'C:\Program Files (x86)\GRASS GIS 7.2.0\grass72.bat'
startcmd = [grass7bin, '--config', 'path']
try:
    p = subprocess.Popen(startcmd, shell=False,
                         stdout=subprocess.PIPE, stderr=subprocess.PIPE)
    out, err = p.communicate()
except OSError as error:
    sys.exit("ERROR: Cannot find GRASS GIS start script"
             " {cmd}: {error}".format(cmd=startcmd[0], error=error))
if p.returncode != 0:
    sys.exit("ERROR: Issues running GRASS GIS start script"
             " {cmd}: {error}"
             .format(cmd=' '.join(startcmd), error=err))
gisbase = out.strip(os.linesep)

os.environ['GISBASE'] = gisbase
grass_pydir = os.path.join(gisbase, "etc", "python")
sys.path.append(grass_pydir)

import grass.script.setup as gsetup

gisdb = r'C:\data\grass'
location = "Country"
mapset = "Topic"

rcfile = gsetup.init(gisbase, gisdb, location, mapset)

import grass.script as gscript
print'grass.script imported'

gscript.run_command('g.gisenv', 'n')

os.remove(rcfile)

EDIT:

After trying to reinstall GRASS GIS several times I decided to download a new GRASS GIS installation file and install that. Now the script runs without throwing the error.

EDIT 2:

I thought the issue was resolved but that is not the case. When I try to use a command as per the tutorial here:

gscript.run_command('r.in.gdal', input = inputVRT, output = outputName, overwrite = True)

It doesn't work and I again get a CalledModuleError:

Traceback (most recent call last):
  File "C:\scripts\10_1_grassGis.py", line 54, in <module>
    overwrite = True)
  File "C:\Program Files (x86)\GRASS GIS 7.2.0\etc\python\grass\script\core.py", line 410, in run_command
    return handle_errors(returncode, returncode, args, kwargs)
  File "C:\Program Files (x86)\GRASS GIS 7.2.0\etc\python\grass\script\core.py", line 329, in handle_errors
    returncode=returncode)
CalledModuleError: Module run None ['r.in.gdal', '--o', '--q', '-o', 'input=C:\\data\\rasters\\bykerneNet_03743.vrt', 'output=bykerne03743'] ended with error
Process ended with non-zero return code -1073741512. See errors in the (error) output.

However, if I try things like:

gscript.gisenv()

It runs fine.

What could be the issue?


texts:
I'm trying to join an excel spreadsheet with one of my layers. I've formatted the spreadsheet to fit QGIS, with basicaly one string and a few integers.
I've saved my spreadsheet in a .csv fromat, and created a text folder saved in .csvt format, with the following code: "String","Integer","Integer","Integer";"Integer","Integer".

I've imported the spreadsheet into QGIS, however, once joined, all my integer become strings, as shown in the image below. Any reason why this happens?


texts:
I would like to create a NDVI Time Series map? How do I do that? The one I probably would like to have similar like this http://proceedings.esri.com/library/userconf/proc02/pap0593/p05931.jpg

I created many NDVI for each or several months from the year of 2002 to 2013?

1) Do I need to add monthly NDVI to create an output monthly NDVI file and then the next step is to take new output monthly NDVI to produce a new output one year file?

(or)

2) Add all the NVDI from one year to generate  a new output for one year file ?

The NDVI were produced from an ASTER.


texts: Converting Degrees, Mintues, Seconds to X, Y to be plotted into ponts

texts:
I need to convert a longs list of Lat/Long (example: N42 30' 45" W122 38' 00") into X,Y data (42.5125, -122.633333) so I can plot it into point shapefile.
Is there a a way to convert multiple records at once from an excel spreadsheet? Maybe an excel formula that I can plug in or a online tool?


texts: How do I create a rectangular feature in ArcGIS Pro given only 2 coordinates and a width?

texts:
I need to create a rectangular feature in ArcGIS Pro that is a representation of a linear tornado path.
I only have the two lat/long coordinates marking the start and points of a line. This straight line is the path the tornado took. I need to 'expand' the line into a rectangle such that this path line and its endpoints are the middle and midpoints of the "short" sides of the tornado path rectangle. I also know the width of the tornado path (how wide the rectangle should be). How do I create an exact rectangle based on this info?
I am using ArcGIS Pro 2.5's editing tools to create the feature in a polygon feature class.


texts:
I've got a polygon feature class with attributes fields
y2010, y2011, y2012, y2013, y2014, country

The 'year' fields contain counts of an event occurrence, which can be 0 or > 0. The 'country' field' is an identifier for country A or B, so the table looks like this:

y2010 | y2011 | y2012 | y2013 | y2014 | country

   0  |  10    |   1   |   0   |   5   |   A
   1  |   0    |   3   |   0   |   5   |   A
   0  |   1    |   1   |  45   |  20   |   B

I'd like to reformat this table to show per year the number of features that fit certain criteria. The criteria are:

More than zero event occurrences for countries A+B;
More than 1 event occurrences for countries A+B;
More than zero event occurrences for country A.

so that the final table using the example above would looks like:

criterium | 2010 | 2011 | 2012 | 2013 | 2014 |
   1           1      2      3      1     3
   2           0      1      1      1     3
   3           1      1      2      0     2

and the table content of columns 2:6 would be the number of features of the original table, that meet the criterium for that specific year.

I've been toying around with the Statistics toolset, but so far without success. Any tips are welcome.

Edit: Opted for Q&D solution by exporting the attribute table to R and wrote a script to calculate the cell values.


texts: Dissolving points produced by GeneratePointsAlongLines_management tool in ArcPy?

texts:
The tool Generate Point Along lines generates points equal to the number of adjacent touching lines. For example, for 3 lines emerging from a point, the number of points generated would be 3.

How do I merge these points to a single point?

I have tried:

arcpy.GeneratePointsAlongLines_management (ORLinks, ORNodes, "DISTANCE", Distance='5000 meters', Include_End_Points='END_POINTS')
arcpy.MakeFeatureLayer_management("ORNodes", "ORNodes")
arcpy.SpatialJoin_analysis ('ORNodes', 'ORNodes', temp3, "","","INTERSECT")

but, this does not seem to work.


texts:
I am trying to compile a list of file extensions that are accepted by GDAL, something like

['.tif', '.png', '.jpg', '.jpeg', ...]

I found a list of GDAL Raster file formats, but the values in the "code" column in that table are not always file extensions.

Is there a similar table that includes all filename extensions?

The goal is to check if GDAL can open a file from a path that is provided by users.


texts: How to draw or mark an area without having any data of that area?

texts:
I'm doing some market research to start a business in Amsterdam. The City consists of 8 city area's which I'd like to map in my map as one of the layers and then add data to it.

I already called to City Hall to see if they had the data readily available for importing into a GIS solution but unfortunately they don't.

The best I could do was find a map with an overview of the city parts and neighbourhoods.

Does anyone have an idea or tip on how to best go about this?
I'm making my map in CartoDB. I'm open to solutions using Mapbox.

I have no programming skills so please bear with me.


texts: Why do imported excel long/lat points not align with basemap?

texts:
Using old plans and maps that I rectified and georeferenced I determined approximate coordinates for certain points. These coordinates are in decimal degrees and are split in two colums latitude and longitude in the relevant excel table. Now I imported the excel table in xls format and set the coordinate system to WGS84 Web Mercator, and while the relative position of the points to each other seems correct, their absolute position is wrong; they are in the middle of the base map and not in Belgium, which does obviously not exist in the ocean near Africa.

What am I doing wrong?

I was under the impression that WGS84 web mercator was the default setting for ArcMap10. I also tried exporting the imported xls into a shapefile. Or am I doing something wrong with the excel since every time I get a 'ArcMap drawing Errors: one or more layers failed to draw. syntax error (missing operator) in query expression ' notification. I tried changing the ',' into '.' in the DD notation but then my computer reads it as numbers instead of coordinates.


texts: Conversion to_string provides NULL in QGIS

texts:
I've been trying to convert a decimal row into a string row using a few methods. Below is my CSV file:

also, note that csv header has been updated too,

State;AverageRainfall;ProbabilityBelow;ProbabilityNormal;ProbabilityAbove

And white space was removed from the headers removed. (including the blank cells)

First method using to string:

The second method using quotes:

Tried the same with integers, a similar result was depicted:

Removing the file name from the column doesn't help:

Using QGIS 3.6.3-Noosa. What I'm doing wrong?


texts: Find index location of a lat/lon point on a raster grid in ArcPy

texts:
This question has been confusing me for a while, and I'd love y'all's help.

I'm trying to find the index location on a raster (in terms of column and row#) of a specified lat/lon point.

My first thought was to find the cell value at the lat/lon point, convert the raster to a numpy array, and find the index of the point (but that's clunky, and requires each cell value to be unique...which isn't true in the raster)

#get values, assuming the extent object is already saved
raster = "raster.shp"
x_coord = extent.XMin
y_coord = extent.YMin
point = arcpy.Point(x_coord, y_coord)
cell_val = arcpy.GetCellValue_management(raster,point,"")
raster_array = numpy.array(raster)
final_index = raster_array.index(cell_val)
print final_index

But that's ungainly, and also requires the cell values to be unique in the raster (which again isn't the case)

Is there any other way to find the cell index of a pixel within a raster that corresponds to a point object?

Thanks!


texts: Save Print/Map QGIS composer view as PNG/PDF using Python (without changing anything in visible layout)?

texts:
I have a QGIS composer/print view opened with all elements adjusted as I wish.
Now I need to print/save/export this as a PNG or PDF file from Python console/Python script.

I do not want to change anything in current layout. Most examples which I found (for example: this), change a map position or size in output PDF comparing to what I see in my current composer view.
I wish to receive exactly the same result as I would get when I click on Print --> Export as Image.

How can I do that? Atlas is not a solution for me.


texts: Arcmap field calculator python script returns no data

texts:
I am running ArcMap 10.2
I have 4 fields with which I am working. many have duplicate entries. in order to find the single correct entry from all the duplicates I am trying to calculate out the erroneous ones. My script runs, but I get no data returned to the GeoloCode field.
Why?

codeblock:

def cal(GeoloCode):
    if ((House_Num >= FROM_ADDR) and (House_Num <= TO_ADDR)) :
        return 88
    else:
        return 1

Expression

cal(!GeolocCode!)


texts: Accessing GDAL from Terminal with OS X Mavericks?

texts:
I successfully installed QGIS on my iMac with OS X 10.9.1 Mavericks using the installers provided by Kyngchaos. However, I'm having troubles accessing GDAL / ogr from terminal. It's giving me the following error message:

  -bash: ogr2ogr: command not found

Does anybody have an idea of how to resolve this?


texts: Adding time value into network model of ArcGIS Network Analyst?

texts:
I am trying to set up a network model in ArcGis 10.3.1 so that I can find a route which will also give the time in minutes. The tutorial guide says this is set  up in the attributes of the network dataset, but within this stage there are various evaluators and I am not sure which one to use. In the attribute table of the feature class I am using I do have a length and time field already calculated.

How to incorporate time?


texts: QGIS Stops Working When I Run the PyQGIS Script the second time

texts:
QGIS Version - 3.0
I developed a script which creates layers and saves it the to the disk. When I run the script, it runs successfully the first time. But when I run it the second time, QGIS crashes. Can someone help me to figure out what could be the problem?

Here is the script.

from qgis.core import *
import qgis.utils

QgsApplication.setPrefixPath("C:/Program Files/QGIS 3.0", True)
qgs = QgsApplication([], False) ###
qgs.initQgis() ###

_layer1 = QgsVectorLayer("C:/AR/Automation Test Script/Radio Signal/MasterLine.shp", "MASTERLINE", "ogr")
_layer2 = QgsVectorLayer("C:/AR/Automation Test Script/Radio Signal/SIGNAL.shp", "SIGNAL", "ogr")
_layer3 = QgsVectorLayer("C:/AR/Automation Test Script/Radio Signal/ClearancePoint.shp", "CLEARANCEPOINT", "ogr")

if not _layer1.isValid():
  print("Layer failed to load!")

else:
    print("Layer was loaded successfully!")

if not _layer2.isValid():
  print("Layer failed to load!")

else:
    print("Layer was loaded successfully!")

if not _layer3.isValid():
  print("Layer failed to load!")

else:
    print("Layer was loaded successfully!")

_writer1 = QgsVectorFileWriter.writeAsVectorFormat(_layer1,"C:/AR/Automation Test Script/Radio Signal/MasterLine12.shp",'utf-8',QgsCoordinateReferenceSystem(),'ESRI Shapefile')
_writer2 = QgsVectorFileWriter.writeAsVectorFormat(_layer2,"C:/AR/Automation Test Script/Radio Signal/SIGNAL.shp",'utf-8',QgsCoordinateReferenceSystem(),'ESRI Shapefile')
_writer3 = QgsVectorFileWriter.writeAsVectorFormat(_layer3,"C:/AR/Automation Test Script/Radio Signal/ClearancePoint.shp",'utf-8',QgsCoordinateReferenceSystem(),'ESRI Shapefile')

qgs.exitQgis()

#_writer = QgsVectorFileWriter.writeAsVectorFormat(_layer,r"C:/AR/RadioSignal/MasterLine123.shp","utf-8",None,"ESRI Shapefile",True)

This is after first run

This is after second run.


texts: Batch reading features from many Geodatabases in R

texts:
I would like to learn to batch import from many geodatabases at once and merge the files together into a new unified shapefile. All of the GDBs and Features have a similar name and file structure. For example GDB1, GDB2, GDB3,...ect. and inside each have shp1, shp2, and shp3.....ect.

From other posts in the group I learned that if these features were not withing GDBs I could simply create a list of the files....

shp_files <- list.files(wd,  pattern = glob2rx("*.shp")

and then use rbind and lapply to merge them together....

merged<-do.call(rbind, lapply(shp_files, rgdal::readOGR))

But how can I alter this code to work if the files are in a GDB or if the bundle of features are in a zipped folder?


texts:
I am using Mapbox's queryRenderedFeatures https://www.mapbox.com/mapbox-gl-js/api/#map#queryrenderedfeatures

I used it successfully to make map symbols interactive
https://s3-us-west-2.amazonaws.com/gomaps.us/clientmaps/CBdevo/index.html

But, I can't seem to get it to work for lines
https://s3-us-west-2.amazonaws.com/gomaps.us/clientmaps/CBMBA/masterplan/index.html

I'm assuming that the code is the same for symbols and lines, but maybe I'm wrong?


texts:
I have GPS data from multiple users and multiple routes and I'd like match them into the actual route using QGIS, because there are some outline points.
My data is a csv file with 2.989 points, containing user-id, route-id, modality of transportation, time and date, x and y coordinates, accuracy, speed and heading from the GPS.
I tried to do that with the Offline-MapMatching plugin, but it seems that it has some bug or something and doesn't work.
When using the plugin, I selected a road layer that I have as a network layer, the GPS points as the trajectory layer and for the trajectory ID I chose the time column from the points. I didn't change the maximum search distance and fast matching type for the plugin to run it.
How can I solve this problem?


texts: How do I make text show up on top of Markers using Carto in TileMill

texts:
I'm attempting to place intersection markers on a trail map. Each intersection has a number. I can't seem to get the text show up above the marker. If I use the opacity, I can see the text is underneath. Is there an option I'm missing?

Here's the output:

And here's the Carto:

#intersections {
  text-name: "[Name]";
  text-placement: point;
  text-fill: black;
  text-face-name: @labels;
  text-size: 12;
  marker-width:10;
  marker-fill:#fff;
  marker-line-color:#333;
  marker-allow-overlap: true;
  marker-placement: point;
  marker-opacity: 0.75;
}


texts:
This code worked fine before but now I'm getting this error:

%%X was unexpected at this time

Here is the code I'm using:

cd /d C:\Users\user\Documents\PhiladelphiaBaseLayers
set ogr2ogrPath="c:\program files\QGIS Essen\bin\ogr2ogr.exe"
set GDAL_DATA=C:\program files\QGIS Essen\share\gdal
for %%X in (*.shp) do %ogr2ogrPath% -skipfailures -clipsrc C:\Users\user\Documents\PhiladelphiaBaseLayers\clipFeature\city_limits.shp C:\Users\user\Documents\PhiladelphiaBaseLayers\clipped\%%X C:\Users\user\Documents\PhiladelphiaBaseLayers\%%X
for %%X in (*.shp) do %ogr2ogrPath% -skipfailures -s_srs EPSG:4326 -t_srs EPSG:3857  C:\Users\user\Documents\PhiladelphiaBaseLayers\clippedAndProjected\%%X C:\Users\user\Documents\PhiladelphiaBaseLayers\clipped\%%X

Went back and forth through the original code but I can't find where I'm off.

I made this code to test the directory for any folders that have "shp" in the name:

import os

files = os.listdir(r"C:\Users\user\Documents\PhiladelphiaBaseLayers")

shpNames = []

for f in files:
    if f[-3:] == "shp":
        shpNames.append(f)
for f in shpNames:
    print f

All files within seem to be legitimate shapefiles

arts_centres.shp
atms.shp
banks.shp
bars.shp
bicycle_parking.shp
bus_stops.shp
cafes.shp
cemetery.shp
fast_food.shp
fire_stations.shp
fuel.shp
hospitals.shp
hotels.shp
landuse.shp
libraries.shp
movies.shp
museums.shp
natural.shp
places.shp
places_of_worship.shp
police_stations.shp
railways.shp
restaurants.shp
roads.shp
schools.shp
stations.shp
swimming_pools.shp
toilets.shp
waterways.shp


texts: Practical way of managing polygons and unions of them

texts:
I have a table in PostGIS with a quite fine partition of a region, which I'd like to dissolve at two higher levels by means of two foreign keys (e.g. municipality and state). The partition is a real one, i.e. the areas don't overlap, they just touch each other.
I thought it'd be nice to do that in views, to follow the "don't insert into the database something which can be calculated from other existing data".
However, a simple

SELECT id_municipality, ST_Union(the_geom)
FROM areas
GROUP BY id_municipality

takes around 22s on a core i7 for "just" around 36000 areas with a POLYGON geometry field.
I indexed the geometry field with a GIST index, but with no luck.

I am now thinking that the only practical way would be to insert a geometry column in the municipalities table and calculate it beforehand.

Unfortunately, I need to share my results to less technical people, so installing experimental packages (like topological extensions) is not practical.


texts:
I would like to multiply all the values within a raster layer with one value, which is written in an *.csv table (it is the only value in this table).

I tried this code (changing Times tool to Minus tool)

>>> import arcpy
>>> from arcpy import env
>>> from arcpy.sa import*
>>> arcpy.env.workspace = "D:/Bachelorthesis"
>>> dgm = ['dgm']
>>> gwf = ['Groundwatertable']
>>> path = 'd:/Bachelorthesis/gis_daten_footprint/Data(1)/GIS_data.mdb/gwf'
>>> value = [row[0] for row in arcpy.da.SearchCursor(path,gwf)]
>>> Gwt = set(value)
>>> print (Gwt)
set([28.124])
>>> outMinus = Minus(dgm,Gwt)

But then this error appears:

 Runtime error
 Traceback (most recent call last):
   File "<string>", line 1, in <module>
   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 3840, in Minus
     in_raster_or_constant2)
   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Utils.py", line 47, in swapper
     result = wrapper(*args, **kwargs)
   File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\sa\Functions.py", line 3837, in wrapper
     return _wrapLocalFunctionRaster(u"Minus_sa", ["Minus", in_raster_or_constant1, in_raster_or_constant2])
TypeError: expected a raster or layer name

But 'dgm' is definitely a raster...
Does anyone know how to make it work?


texts: FME Keeping attributes of individual polylines joined together by LineJoiner transformer

texts:
Quick Q: I've seen some other questions around this but not quite what I'm looking to do. Hopefully it's simple.

I have a series of lines that I want to join together using the LineJoiner transformer. So far so good and this works nicely...BUT...each polyline has an ID unique to that line segment.

I would like to use the LineJoiner to chain polylines together but I want to keep the combined attributes (ID) as a list associated with the elements (either comma seperated or as a seperate list). The reason is that at the end I'd like to be able to see a polyline and know it's been made up of ID's 2,3,4.

Is this possible in FME or do I need to look elsewhere?


texts: Incorrect transformation when using pyproj

texts:
I am using pyproj to work with shp files and am relatively new to GIS. I am tryign to visualize the shp data, but am having a problem as I am incorrectly transforming the coordionates.

Here is an example from my code:

import shapefile
shp = shapefile.Reader(filePath)  # .shp file
rec = shp.shapeRecords()[0]
x, y = rec.shape.points[0]

Now, the resulting x and y are in utm format. For plotting, I want to convert to (lat,lon). I am using the following:

import pyproj
p = pyproj.Proj(init='EPSG:3857', proj='utm', zone=36N, ellps='WGS84')
p(x,y,inverse=True)

But this gives me incorrect (lat,lon) coordinates. the point in question is close to the Syrian-Israeli border, and has the raw value of (3970162.65625, 3853896.40625). But running p(x,y,inverse=True) gives the (lat,lon) coordinates of 67.809, 29.709, which turns out to be around Finland.

Any ideas why this transformation is wrong? I checked for the zone in here. And I am pretty sure the other parameters are correct as well. The shp file was created by exporting from ArcGIS desktop's ArcMap.

Edit in response to comment:

Here is the .prj file:

PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",
GEOGCS["GCS_WGS_1984",
DATUM["D_WGS_1984",
SPHEROID["WGS_1984",
6378137.0,298.257223563]],
PRIMEM["Greenwich",0.0],
UNIT["Degree",0.0174532925199433]],
PROJECTION["Mercator_Auxiliary_Sphere"],
PARAMETER["False_Easting",0.0],
PARAMETER["False_Northing",0.0],
PARAMETER["Central_Meridian",0.0],
PARAMETER["Standard_Parallel_1",0.0],
PARAMETER["Auxiliary_Sphere_Type",0.0],
UNIT["Meter",1.0]]


texts: Polygonize vegetation area from Landsat7 RGB using QGis / Bilko?

texts:
I'm trying to extract the areas with vegetation as a polygon-shape from a landsat7-image, 30m resolution, RGB, using QGis 1.8 or Bilko or (hopefully in the near future) ArcGIS. The area is in Northern Peru, dry season, so basically I want to extract all the shades of green.

I have never really worked with raster data and don't know what to do. I have found out there are "vegetation indices" which can be calculated, and have found Bilko as a software I could do this with. However, if I got it right I need more bands than RGB to do this - is this true? I'm not sure.

On this website I have also found threads on how to calculate vegetation area in QGis - using GDAL tools or SEXTANTE: Calculating area of vegetation in raster file?

However I would like to not only calculate but extract the data as a polygon.

What would be the best way to do this? I've never worked with Bilko, GDAL or SEXTANTE, of course I will go into that and try to learn, but I would like to know what will work before I spend hours and hours only to find out I can't do it...

P.S.: I will probably have ARCGis in a while, so if it is much easier with ARCGis I might wait for that...


texts:
I wrote a Parser for a GML file from the OS Mastermap with geotools 8.4. With one part I still have a problem. This sniplet of the xsd:

<complexType name="RoadLinkType">
    <complexContent>
        <extension base="osgb:AbstractFeatureType">
            <sequence>
                <element name="polyline" type="gml:GeometryPropertyType"/>
                <element name="directedNode" type="osgb:directedNodeAssociationType" minOccurs="2" maxOccurs="2"/>
            </sequence>
        </extension>
    </complexContent>
</complexType>

The part from the gml file:

<osgb:RoadLink fid='osgb4000000023204016'>
    <osgb:polyline>
        <gml:LineString srsName='osgb:BNG'>
            <gml:coordinates>516609.000,257678.000 516615.000,257733.000 516618.000,257786.000</gml:coordinates>
        </gml:LineString>
    </osgb:polyline>
    <osgb:directedNode orientation='-' xlink:href='#osgb4000000027916595'/>
    <osgb:directedNode orientation='+' gradeSeparation='1' xlink:href='#osgb4000000028203009'/>
</osgb:RoadLink>

It gets parsed nearly right, i.e. the polyline is correct (I get a class com.vividsolutions.jts.geom.LineString) but I only get one (instead of two) directedNode.

My code:

GML gml = new GML(Version.GML3);
CRSAuthorityFactory crsFac = ReferencingFactoryFinder.getCRSAuthorityFactory("EPSG", null);
CoordinateReferenceSystem osgbCrs = crsFac.createCoordinateReferenceSystem("EPSG:27700");
gml.setCoordinateReferenceSystem(osgbCrs);
SimpleFeatureIterator iter = gml.decodeFeatureIterator(in);
while (iter.hasNext())
{
    SimpleFeature feature = iter.next();
    System.out.println("id: " + feature.getID());
    System.out.println("polyline: " + feature.getAttribute("polyline"));
    Object dnode = feature.getAttribute("directedNode");
    System.out.println("directedNode: '" + dnode + "'\t" + dnode.getClass());
}

And the output:

id: osgb4000000023204016
polyline: LINESTRING (516609 257678, 516615 257733, 516618 257786)
directedNode: '{orientation=-, href=#osgb4000000027916595}' class java.util.HashMap

What am I doing wrong? Why don't I get both directedNode? How can I do the binding to real Java Objects?


texts: Count duplicates in a field for sequential unique ID

texts:
I am trying to create a unique ID for my street feature class. (My unique ID starts with first 3 letters of the Street Name, 2 letters of the Street type, and 3 letters of the City Name Abbreviated.)
Where I am stuck is that I am trying to get a number to increment on the segments that would be duplicated.  For example, Glenacre Drive has three segments. The UniqueID would be GLEDRHQF01, GLEDRHQF02, GELDRHQF03.
I cannot figure out how to create sequential numbers not based off of the ObjectID sorting... For example, I used the following Field Calculator but realized it is looking at the ObjectID sort and my streets aren't entered in order...
Pre-Logic Script Code:
prevFieldValue = ''
counter = 1
def GetDuplicateCounter(myFieldValue):
  global prevFieldValue
  global counter
  if myFieldValue == prevFieldValue:
    counter += 1
  else:
    counter = 1
  prevFieldValue = myFieldValue
  return counter

=
GetDuplicateCounter( !UNIQUEID! )


texts: Using ModelBuilder to iterate through and export to single layer

texts:
I have a geodatabase that has multiple layers, the layers all have some of same attribute names. I am trying to use ModelBuilder to iterate through the GDB and export all polygon features that have the same attribute name into a single layer.
For example, I have a buildings, parks, and structure layers in a GDB all three have the same attribute names "facilityNumber", "facilityIDNumber", and "facilityName" among others that are different for each layer. I want to export just the shape and the attribute information that is the same across all features into one single polygon layer.


texts: How to create an accurate Tissot Indicatrix?

texts:
A Tissot Indicatrix is useful method for communicating at a glance the kinds of distortion a given projection is prone to (in the figure below, each of the red circles occupies the same area). I've been told that popular methods for generating TI's have their own problems, to the point of being sometimes woefully inaccurate.

What is the problem with the popular methods, and what is most correct way of generating a TI that is accessible to your average GIS dude(ette)?


texts: Changing longitude,latitude to latitude,longitude in OpenLayers 3?

texts:
Openlayers is representing the coordenates as longitude,latitude.

Can I change the default format longitude,latitude to latitude,longitude?

By example, this will not work, because I have this array as [lat, lng] and LineString is taking it as [lng,lat]:

a = [
    [-21.929054, 64.127985],
    [-21.912918, 64.134726]
    ];
var lineString = new ol.geom.LineString(a);


texts: DBMS table not found error adding PostGIS layer in ArcMap?

texts:
I imported some shapefiles to a postgis database. In QGIS it works fine, however in ArcMap 10.4 not yet. I added postgreSQL as a database connection in ArcMap. I can also see the postgis layers in the ArcCatalog. But whenever I try to load them to the ArcMap document, I get the error:

  DBMS table not found


texts: Changing scale of one map sheet within map series in MPS Atlas?

texts:
I am using MPS Atlas and have a large map series (almost 300 map sheets) with the Areas Of Interest defined by the polygons of a feature class with a 100 m buffer on each polygon. This works well for most of the map sheets. However, for 7 of the map sheets the scale needs to be larger (i.e. they need to be zoomed out further than the extent of the polygon).

How can I change this with out changing the extent for the rest of the map sheets?

I have tried making the map sheets I want to change "single" (i.e. not shared) but when I do this I cannot find an option to change or reset the scale.

Does anyone know how to do this?


texts: Legend item resolution issues in ArcGIS Javascript app

texts:
When I create a pdf map using my JavaScript API app, the legend items have messed up resolutions when there are more than a few swatches in the legend layers:

On the left: image of print with poor resolution on swatches, on the right: image of .mxd template
 .

The only way to correct this is seems to be to make the legend one column at the service level:

On the left: image of print with good resolution on swatches, on the right: image of .mxd template
 .

The only parameter in the 'Web Map as JSON' string that seems relevant is the DPI, but changing that has no effect. I have also messed with the Legend properties dialog box in the corresponding mxd to no avail. The bug seems to occur when the legend adjusts the number of columns. This issue is more than a year old for me, with no luck googling the heck out of it... Any ideas?


texts: Polygon made from raster extent does not display exactly same as raster in Google Earth Pro?

texts:
Polygon made from a raster extent does not display exactly the same as the raster in google earth pro.

Does anyone else have this problem?

 new_crs<-crs(dem_10)#get latlon
 r2<-projectRaster(r1, crs=new_crs)# project to latlon
 p <- as((extent(r2)), 'SpatialPolygons')# get extent polygon
 KML(r2, file="example_data.kml")# write raster to kml
 KML(p, file="example_data_extent.kml")#write polyogn to kml


texts:
My purpose is to use the Mxd2Qgs plugin in order to be able to use mxd files with QGis.

I have downloaded the Mxd2Qgs file and unzipped it. I put it in the folder:

C:\Program Files\QGis Lyon\apps\qgis\phyton\plugins

However, when I open the Plugin Manager in QGis 2.12, the Mxd2Qgs plugin does not appear. How do I resolve this? I have only very basic informatic knowledge (I read similar questions but it did not help me resolve the problem).


texts: How to add a url(having a webseal) in wfs service of QGIS?

texts:
I am using QGIS 2.18 version.I am trying to add a new connection in wfs. That url is having a webseal. While I am trying to add a new connection in the WFS, the log message i showing like this:

Download of capabilities failed: Host requires authentication

I have given the username and add basic authentication also. Is it Possible with other versions?


texts: I want to assign class by thematic layer in ecognition. I have digitized road edges and I want to classify the roads in my 'road class' using the shapefile (polyline) of my roads. How do I do this?


texts: Automatically detecting relief in flat areas from hillshape for marine bathymetry using ArcGIS Desktop?

texts:
I am trying to identify areas of interest in the ocean around 30-150 m that have some "topographic relief" compared to their surrounding areas.  I did this manually by looking at the hillshade and circling areas that look like they "bulge out" a little big.  However, I know this is not that advanced and was wondering if there was an automated method for detecting relief and creating new polygons from this detection.

I am a beginner.


texts:
I need to download the osm2pgsql for Windows

i tried from the path which is given in wiki

http://customdebug.com/osm/osm2pgsql.zip

and also from below link

http://tile.openstreetmap.org/osm2pgsql.zip

Both are not working

Please share me if any one has this for Windows 64Bit


texts: Edit and Update Existing Features in a Collector App using Survey123 Forms

texts:
I have an existing polyline dataset that needs to be reviewed and have certain attributes populated by field personnel. We have a hosted feature layer in Collector, and are attempting to use a Survey123 form to do the review as the form has much greater customization and further capabilities than simply editing in Collector.

The ideal workflow would be:

Click feature in Collector to bring up a custom pop-up
Click custom url in pop-up to launch Survey123 field app with certain fields populated.
Submit survey to populate/overwrite attributes in the existing feature.

url: arcgis-survey123://?itemID=123455abcdefg&mode=edit&globalId={GlobalID}&objectId={OBJECTID}&field:PermId={PermId) etc...

Where I am getting hung up is apparently you cannot use the mode=edit parameter with the Survey123 Field App, only a Survey123 web form? Trying to use these parameters in the Field App url still is adding a new record in the underlying table rather than updating it. Is there a way to either have the form submission overwrite an existing record or open survey123 straight in the inbox rather than a new survey?

We want to use the map capabilities and interface of the collector app with the form customization and data validation abilities of the Survey123 app. We cannot use a web form to edit the features (frustratingly, since that seems to work fine) as some of the personnel will be working in areas without data connection. We need the offline capabilities of the Field App.


texts: I am an user of QGis 1.9.0. When I go to Fetch plugin option Sextante plugin is broken seems what is the problem?Even I agin downloaded the plugin in Qgis official resportary but it displays me the like below:


texts: AIS Data Gaps - plotting as lines in QGIS

texts:
I am currently attempting to plot AIS data from our on board AIS receiver in QGIS. (https://en.wikipedia.org/wiki/Automatic_identification_system)

My current processing steps are :-

Import the AIS data as a CSV file containing Time, Vessel Name, MMSI No, Lat & Lon.
Save imported data as a shapefile
Remove AIS positions outside area of interest using select feature by polygon, invert selection and delete selected.
Convert points to paths using order by time and group by MMSI No

This produces paths of the vessel track.

My problem is that vessels occasionally have a large gap in their AIS data that leads to un-realistic paths, i.e passing over land.

Is there are way to sort by MMSI number (Unique to each vessel) and then assign a unique cluster number to the AIS record, based on the distance between the AIS record positions for each MMSI?

i.e  If a vessel has a gap in it's MMSI data of greater than 5000m assign a new cluster number to that record.

I would need to do this on a vessel by vessel basis, so that only the the AIS data for each specific vessel is examined.

If I then use points to paths and group by cluster number, it should then draw the path for all AIS data points with a gap of less than 5000m.

I followed  @Spatialthoughts’ instructions and they worked perfectly.

Removing segments over 5 km long has improved the readability of the information and produces a more realistic plot of the AIS data.


texts: Creating geotagged twitter maps in cartodb for one country instead of entire world

texts: I know we can create a geotagged twitter maps using cartodb http://docs.cartodb.com/tutorials/connect_twitter_datasource/#twitter-api-data-options. But by default, the cartodb will create a twitter map for entire world. Is there any way to restrict for one particular coutry?


texts: Applying multiple CQL_FILTER to single GeoServer LayerGroup?

texts:
I've created a page with OpenLayers 3 that adds a WMS Tile layer to the map served up by a local GeoServer (v2.6) instance.  As it stands adding layer by layer things work just fine:

var dynamicMap = new ol.Map({
  target: 'myMap',
  layers: [
    new ol.layer.Tile({
      source: new ol.source.TileWMS({
        url: 'http://demo.opengeo.org/geoserver/wms',
        params: {
          'LAYERS': 'ne:NE1_HR_LC_SR_W_DR'
        }
      })
    }),
    new ol.layer.Tile({
      source: new ol.source.TileWMS(({
        url: 'http://localgeoserver:8080/geoserver/wms',
        params: {'CQL_FILTER':"[(USAGE IN ('H','B')]", 'TILED': true, 'layers':'WORKSPACENAME:LAYER1','WIDTH':'256','HEIGHT':'256' },
        serverType: 'geoserver'
      }))
    }),
    new ol.layer.Tile({
      source: new ol.source.TileWMS(({
        url: 'http://localgeoserver:8080/geoserver/wms',
        params: {'CQL_FILTER':"[SYMBOL IN ('A','B','C')]", 'TILED': true, 'layers':'WORKSPACENAME:LAYER2','WIDTH':'256','HEIGHT':'256' },
        serverType: 'geoserver'
      }))
    })
  ]
});

I plan to add several more layers to GeoServer, and rather than making WMS requests for each layer, I want to put them all in a Layergroup in GeoServer.  I did that --

new ol.layer.Tile({
      source: new ol.source.TileWMS(({
        url: 'http://localgeoserver:8080/geoserver/wms',
        params: {'TILED': true, 'layers':'WORKSPACENAME:LAYERGROUP1','WIDTH':'256','HEIGHT':'256' },
        serverType: 'geoserver'
      }))
    })

and again all the layers return as I expect them to, with one exception - I'm not sure how to apply the necessary CQL_FILTER to my layergroup.  When I added the first CQL_FILTER [(USAGE IN ('H','B')] when adding the LAYERGROUP1 layergroup as a layer, it failed because USAGE was not found in the LAYER2 layer.

I tried examples I saw where people said that each CQL_FILTER needs to be delimited by a semicolon, but that didn't seem to work either.

Is my end goal even possible going this route?

Is there a way to apply CQL_FILTERS when creating my layers in GeoServer?

Each layer is simply referencing an Oracle table, so should I go so far as to just create each layer with SQL and just apply my filters in the WHERE clause there?


texts: Using MapBasic, is there a method to compress GeoTIFF after programatically creating the GeoTIFF?

texts:
I need to compress the file sizes of GeoTIFFs that are created when running the MapBasic application. In a .Net application, I transform the GeoTIFF into a compressed JPG using ImageMagick; which works perfectly. However, the original GeoTIFFs are required for sharing and download. Adding to the difficulty is that the existing resolution and dimensions must be maintained. I am trying to perform the compression within MapBasic so that an additional application just to compress the GeoTIFFs will not need to be created.
Here is a sample of the current code:
Save Window windowId As reportPath & imgFileName & ".tif" Type "GEOTIFF" Width 7.49826 Units "in" Height 9.5816 Unites "in" Resolution 1152


texts: QGIS: Getting the right colour symbology legend colors because of a single color by expression

texts:
This question has been related to this one:
Making colours from table with colour percentages/RGB-codes
In that scenario I have made sure that my layer has the right colours based on a expression. This is just because this certain dataset has a "standard" we all have to follow.
The thing is this is based on a single symbol expression (aka the hexcode rgb). The colours are right, however when I want to construct a legend by the "ouderdom" layer, the colours of the legend are wrong. This is of course because when I try to categorize it by that field, in the legend random colours are chosen, causing the legend to display the wrong colours.
So to clarify my symbology of the layers needs to based on the "ouderdom" field categorization but the colouring needs to still happen the same expression (so no random colours), how do I do this?
To illustrate what I mean I will ad some screenshots here:
Now this is the setting I use to set the right colour to my features in the layer:

In this second screenshot I will get a clear coloring like this:

Which matches the standard map of the dataset seen here, this is the standard of the legend everyone has to follow in terms of coloring:
https://www.dov.vlaanderen.be/sites/default/files/pfiles_files/tertiairkaartA3_0.pdf
The problem is, as you may have noticed, the legend highlighted in this screenshot here:

So now I hope you see and understand what the problem is.


texts: In PostGIS why does spatial index still matter for a table with only points?

texts:
According to this Boundless link, it turns out that in PostGIS, spatial query on table with spatial index uses two pass system: with the first pass approximates spatial relationship with bounding box.

Now i was wondering if such approach will be effective for a PostGIS table with geometry type "ST_POINT", as bounding box of a point IS essentially the point itself.

With this in mind, I run an ST_INTERSECT query by Jmeter to 2 tables with exactly the same set of points, and exactly the same intersect BOX. The only difference is one table has spatial index and the other don't.

The result surprises me: the one with spatial index performed much faster.

In my machine, the average response times are:

with spatial index: 911 ms
without spatial index: 62038 ms

Now the question: why does spatial index still matter for a table with only points?


texts: How to Obtain the feature attributes from feature collection

texts:
I am able to obtain the features of the layer which is selected by the user in the Layers table and then intersected by a user defined polygon, using e.features. I was able to obtain the attribute values of the feature as a property, like this e.features[0].attributes.Family or e.features[0].attributes.Species. The e.features[0].attributes object looks like this, Object { Family="Pinaceae", Species="Pinus cembroides"}.

Now the problem that I am facing is when a user selects a different layer, say states, the attributes properties change, the object now looks like this,
Object { STATE_NAME="Arizona", STATE_FIPS="04", SUB_REGION="Mtn", more...}. So now I need to fetch the attribute data by again changing the properties, like e.features[0].attributes.STATE_NAME.

So I am trying to create a common function where I can obtain the properties of the attributes object and then directly fetch the values for any layer. This is the code that I am trying to implement,

var attrData = [];
for (var i = 0; i < e.features.length; i++) {
    //Work on how to fetch properties for all the layers attributes
    feat = e.features[i];
    for (prop in feat.attributes) { alert(prop); } //this gives all the property values
    var values = "";
    for (prop in feat.attributes) {
      value = (value === "") ? feat.attributes.prop : value += feat.attributes.prop; //But feat.attributes.prop is undefined

    }
  attrData.push([value])
}

As you can see, for (prop in feat.attributes) { alert(prop); } gives me all the property names but when I assign this prop to feat.attributes.prop to fetch the value, it returns undefined.

Not sure where am I going wrong? This is my application, http://128.196.142.12/geo/test/test_new.html. Click on the 'draw poly' in the toolbar and then select 'Pima Pine' layer in the Map Layers column and create a polygon and intersect it with the blue stars (Pima Pine), you should get the table with the info of the stars which intersect with the vector polygon. Now I am trying to replicate the same for the USA States and that is where I am failing.

How should I change the code to fetch the attribute values for all the layers?


texts: Duplicating Features based on dates from CSV file

texts:
Recently I was able to get some data that has some time elements attached to it. This data is for estimated house values by census tract and replicates itself every month when available. The purpose of this would be do some analysis on housing prices over time and create a time lapse map.

After doing a bit of research it seems like these are two different tables. The time lapse table would need to have separate features for every z_date so it would have to be a many-to-many join? My first attempt at doing this was trying to do a traditional join by right clicking the census tract shapefile and joining it to the csv file that contained all of the data. This did not work in my favor as it only returned the earliest dates from the csv file the correlated with the census tract FIPS code. When I go to join it it comes back with only the earliest record for all of the census tracts.

The analysis file would have to be a file that had values across each row for a given period of time (average value by quarter) and then create a percent change by acre column for each year. Is there a better way to do this than grouping in SQL using

  WHERE z_date between 01/01/1996 and 12/31/1996

or

  WHERE z_date like '%1996'


texts: Problems using the QGIS Db Manager and accessing a PostGIS file with a JSON field content to plot up JSON attribute content

texts:
Using the QGIS (Madiera) Db Manager window and accessing a PostGIS file to try and plot up some of the JSON file content on to an existing QGIS basemap.  I am not a QGIS (or GIS) superuser!
I have used the SQL query in the DB Manager select coordinates, json_record ->> 'monitor_location' as location from field_observations;
to extract from a JSON field "field_observations".  This SQL query works, revealing the relevant attribute in the Db Manager viewer.
I then tried to use the “Load as new layer” feature to load the queried dataset, THIS NOW WORKS though there is no CRS assigned.  The preview works, though with the data appearing in an odd-ball place outside my NAD10UTM 40 project map space.
Any suggestions on what to do next to get a visible map?


texts: QGIS crashing when using 'save as' to save a new projection

texts:
Every time I try to save a layer with a new projection using 'save as' my QGIS crashes. The new files are there but they seem false. Please see two screenshots attached, one before and one after I used 'save as'. I also included a screenshot of my 'save as' page. My new projection is supposed to be centered around lat 70° but it does not seem that way. In Addition about 10 rows of data (in the form of rectangles in the screenshot and rows in my attribute table) are missing.

Please note this is also a screenshot of the areas when i changed the project projection to my custom projection, the areas still look off:

I seem to be loosing data when clipping my shape file at 89.9°.

To do this I created a new vector layer to use as a clip layer with an upper limit at 89.9°. I imported it and used the function: Vector --> geoprocessing tools --> clip. Here is a screenshot of what I entered into this function. Any ideas what I could be doing wrong? P.s. There is no 'processing completed' message at the end of this it appears to get stuck in the clipping process.


texts: QGIS: Waypoint name (symbol label) in two rows

texts:
I have couple of longer waypoints names that I would like to put in two rows with center formatting but I just don't know how? I tried to put  tag but it doesn't help...

For example I have the waypoint with name 'one two three five six seven' and I want to make it like this:

one two
three five six seven


texts: 6 degree of freedom helmert transformation

texts:
I'm looking to do datum conversions for points with six degrees of freedom, that is, they are described by a heading, pitch, and roll in addition to their position information.

I know how to convert the orientation vector from a local NED or geodetic frame to the ECEF frame (using the method in this paper), but I don't know how to handle the orientation change when I use helmert transforms to convert between the different datums in ECEF space (i.e. WGS84 to NAD83).

So my question is, how can I perform a 6-DOF helmert transform, which preserves orientation information across the datum transformation? Is orientation affected by the transformation at all?


texts: Use of Google Earth Pro derived area and perimeter for a thesis?

texts: I have ground-truthing data of centroids for 23 sites and their surrounding land uses, but no money to purchase remotely sensed images of this remote region of Latina America. I would use the polygon tool option from Google Earth Pro to trace the land uses (buffer 100m, 200m, 500m from centroid) and copy the estimated area and perimeter given by the Google Pro polygon tool. I have searched extensively for free images for three years and have not been able to find any. Would this infringe Google terms of use? These data will be used in a thesis.


texts: How to work with Multivalues in arcgis script?

texts:
I am  really beginner in arcgis scripting. I was unable to find a solution of my problem, so I'm writing here.

My problem:

On my input I want some shapefiles. When I load them from ArcGis Tool, I obtain something like this:

'C:\Users\Steve\Disk Google\BP-Kaiser\In\Cesta_Clip.shp';
'C:\Users\Steve\Disk Google\BP-Kaiser\In\Pesina_Clip.shp';
'C:\Users\Steve\Disk Google\BP-Kaiser\In\SilniceDalnice_Clip.shp'
'C:/Users/Steve/Disk Google/BP-Kaiser/In/Cesta_Clip.shp'

I want to have these shapefiles in a list, so I convert it using split() function to look like this:

list = ["'C:/Users/Steve/Disk Google/BP-Kaiser/In/Cesta_Clip.shp'", "'C:/Users/Steve/Disk Google/BP-Kaiser/In/Pesina_Clip.shp'", "'C:/Users/Steve/Disk Google/BP-Kaiser/In/SilniceDalnice_Clip.shp'"]

And then I would like to work with it.

But if I use just a simple function:

arcpy.CopyFeatures_management(list[0], outWorkspace+"/copy.shp")

it doesn't work. But if I enter the path to the file directly, it works.

I think that it has to be something simple, but like I said I am really beginner.

If you wold like to help me, but you will need more information, write me please.


texts: ArcGIS: Least-Cost Path between points; Source Destinationissue

texts: I am trying to find the least-cost path between 2 points. The problem I am having is that if I use point A as source and point B as destination, I get a value, but when I do the opposite i.e., use point B as source and point A as destination I get another value for the least-cost path. I need to have the same distance whatever source or destination I am using; is that possible? Is there a way to ask ArcGIS to find the least-cost path that would result in the same distance whatever the source and destination as long as I am using the same two points? thanks


texts: Doing expression for project relation in PyQGIS

texts:
I've build an relation between two vectorlayers. Now i want to write an expression for the relation. The standard expression which is given, is to filter features from the first field. Did anybody have an idea how to change this with python?

Code example:

rel = QgsRelation()
rel.setReferencingLayer(str(map1.id()))
rel.setReferencedLayer(str(map2.id()))
rel.addFieldPair('map1field','map2field')
rel.setId('relation')
rel.setName('to_do_list')
QgsProject.instance().relationManager().addRelation(rel)


texts:
I am working on creating a OSM Server and want to display the correct map of India with disputed boundaries. I have followed the following link: http://xsce.org/wiki/generating_map_tiles#disputed_borders to correct India map using geojson file. But I am getting below error messages in my rendering logs. Can anyone help in resolving this issue?

renderd[23400]: An error occurred while loading the map layer 'ajt': OGR Plugin: cannot find layer by name 'OGRGeoJSON' in dataset '/usr/local/share/maps/style/osm-bright-master/geojson/boundary_claimed_by_india/map.geojson'  encountered during parsing of layer 'india_boundary' in Layer at line 1041 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'

renderd[23400]: An error occurred while loading the map layer 'ajt': OGR Plugin: cannot find layer by name 'OGRGeoJSON' in dataset '/usr/local/share/maps/style/osm-bright-master/geojson/boundary_claimed_by_india/map.geojson'  encountered during parsing of layer 'india_boundary' in Layer at line 1041 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'

renderd[23400]: An error occurred while loading the map layer 'ajt': OGR Plugin: cannot find layer by name 'OGRGeoJSON' in dataset '/usr/local/share/maps/style/osm-bright-master/geojson/boundary_claimed_by_india/map.geojson'  encountered during parsing of layer 'india_boundary' in Layer at line 1041 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'

renderd[23400]: An error occurred while loading the map layer 'ajt': OGR Plugin: cannot find layer by name 'OGRGeoJSON' in dataset '/usr/local/share/maps/style/osm-bright-master/geojson/boundary_claimed_by_india/map.geojson'  encountered during parsing of layer 'india_boundary' in Layer at line 1041 of '/usr/local/share/maps/style/OSMBright/OSMBright.xml'


texts:
My job requires that I produce several print maps every day, for a variety of different users. It is usually pretty straight-forward, but sometimes I have to make a map that uses rather confusing/overlapping data for people who are unfamiliar with the subject matter.

Is there a method to have the colour of a label automatically match the colour of the polygon indicated in the symbology tab? I realize I can select 'Define classes of features and label each class differently' and 'Get Symbol Classes' under the Labels tab, but I still have to manually enter the CMYK values for each symbol class, as each label is automatically set to black.

If there was a method to set the colour automatically, my job would be much less tedious.


texts: Which files format should I favour for import into PostgreSQL+PostGIS?

texts:
Although not a newbie, I am not a GIS expert, so I can use your help.

I have this GIS project, based on a PostgreSQL+PostGIS database.
I import GIS data from many sources. Most of them offer several formats and ways to retrieve the data I will be needing, also because the data is often available from several servers (it's mostly governmental open data).

So far, we had been usually downloading manually shapefiles and importing them into the PSQL database using the shp2pgsql command line tool. We should also use some raster data in the future.

I intend to improve this by automating the whole procedure using Python scripts.
That led me to find out more about WFS servers. Most of the sources from which we need to retrieve data have WFS servers; it is often also possible to download files (.csv or zipped shapefiles) from the web sites. Retrieval from these sources can be automated either way (WFS or from the web sites).

I understood that when doing the http request to the WFS servers, there is the possibility to specify the wanted file format using the parameter outputFormat, although the queried server may only support certain formats.

Question 1: I have not found which file formats may be requested and how to interrogate the server about which file formats it supports.

Question 2: assuming that I have the choice, after stumbling onto the web site http://switchfromshapefile.org/, I know wonder which format I should favour for this use (import into my database). I understand that:

Obviously a PSQL database dump is the best, but seldom available.
Shapefiles not a great choice due to truncated names, absence of CRS definition, restricted number of supported datatypes, unspecified character set (we have a lot of é,è,ç,à here), etc...

But then, that still leaves me to choose between .gml, .gpkg, .kml, .geojson, ...

From what I read, .gpkg and .gml seems the best choice if I want to avoid losing information, albeit with the cost of increased file size.
What would you use?


texts: Choosing feature classes by wildcard from file geodatabase in ModelBuilder?

texts:
This must be easy and I have done it a few years ago (I can't recall if it was using a custom python script inside the model) but is there a way to select certain feature classes from a geodatabase to then feed it to another tool (merge in this case)?
I tried Model Only tools - Select Data and Iterators - Feature Classes and Datasets to collect values but it doesn't work...

So I need to select all the files that have r_End_ for one merge and r_Start_ for a seperate merge.
I had a look at Extracting specific feature class from dataset in ModelBuilder? and https://geonet.esri.com/thread/102308 as well as the arcmap help.


texts: Adding feature object to PostGIS feature type from GeoServer

texts:
I am trying to set up a system that will post feature objects periodically to a PostGIS db and display them via GeoServer and OpenLayers on a web site.
Is it possible to add the new feature objects to the db via GeoServer using the REST API or geoserver-manager (my system is Java based) or will I have to add the features directly into the PostGIS db?


texts: Creating installer for ArcGIS toolbar with Visual Studio 2005?

texts:
I'm working on a small toolbar for ArcGIS 9.2, Windows XP, using Visual Studio .NET 2005. I followed Esri's instructions for creating an installer: http://edndoc.esri.com/arcobjects/9.2/NET/0df20605-b457-42d6-b63c-341a3824474a.htm. I can install the resulting .msi or .exe file on the workstation where the toolbar was developed, without being an admin. When I or our admin. tried to install it on other computers in the office, the install failed with an error code 2869. A web search ('install error 2869') returned mostly vague, unhelpful, and irrelevant results.

Has anyone else experienced this error when creating an installer, and how did you resolve it?


texts:
https://jsfiddle.net/324h2d9q/87/

I tried to make a map with Leaflet where I have multiple layers where I sort things via category, but it seems that it is not working... If I click on ,,Muzee", it should run the addMuzeu1() function, but it doesn't work...any suggestions?

function addScoala1() {
   for(var i=0; i<scoala.length; i++) {
    var greenIcon = new L.Icon({
  iconUrl: 'https://cdn.rawgit.com/pointhi/leaflet-color-markers/master/img/marker-icon-2x-green.png',
  shadowUrl: 'https://cdnjs.cloudflare.com/ajax/libs/leaflet/0.7.7/images/marker-shadow.png',
  iconSize: [25, 41],
  iconAnchor: [12, 41],
  popupAnchor: [1, -34],
  shadowSize: [41, 41]
});
    var marker = L.marker( [scoala[i]['latitudine'], scoala[i]['longitudine']], {icon: greenIcon}).addTo(map);
    marker.bindPopup( "" + scoala[i]['scoala']+"
Detalii:" + scoala[i]['detalii'] + "
Telefon: " + scoala[i]['telefon']);
   var scoala = JSON.parse( '<?php echo json_encode($scoala) ?>' );
   }
  }
   function stringToGeoPoints( geo ) {
   var linesPin = geo.split(",");

   var linesLat = new Array();
   var linesLng = new Array();

   for(i=0; i < linesPin.length; i++) {
    if(i % 2) {
     linesLat.push(linesPin[i]);
    }else{
     linesLng.push(linesPin[i]);
    }
   }

   var latLngLine = new Array();

   for(i=0; i<linesLng.length;i++) {
    latLngLine.push( L.latLng( linesLat[i], linesLng[i]));
   }

   return latLngLine;
  }

  function addZona1() {
   for(var i=0; i < zona.length; i++) {
       //console.log(zona[i]['geolocatii']);
    var polygon = L.polygon( stringToGeoPoints(zona[i]['geolocatii']), { color: 'red'}).addTo(map);
    polygon.bindPopup( "" + zona[i]['nume']);
   }
   var zona = JSON.parse( '<?php echo json_encode($zona) ?>' );
  }

function addMuzeu1() {
   for(var i=0; i<muzeu.length; i++) {
    var marker = L.marker( [muzeu[i]['latitudine'], muzeu[i]['longitudine']]).addTo(map);
    marker.bindPopup( "" + muzeu[i]['muzeu']+"
Detalii:" + muzeu[i]['detalii'] + "
Telefon: " + muzeu[i]['telefon']);
   }
   var muzeu = JSON.parse( '<?php echo json_encode($muzeu) ?>' );
  }
var markersA = [];
var markersB = [];
var markersC = [];
var groupA = L.layerGroup(addScoala1);
var groupB = L.layerGroup(addMuzeu1);
var groupC = L.layerGroup(addZona1);
//background tile set
var tileLayer = {'Harta' : L.tileLayer('https://api.mapbox.com/styles/v1/mapbox/streets-v10/tiles/{z}/{x}/{y}?access_token=pk.eyJ1Ijoicm9iaTEwMTAxMjk4IiwiYSI6ImNqNjZocndibTE2NmQzM3F4YWRzMDBubWkifQ.ccMsQCNcMHOG4ucR6ogCcg', {
            attribution: 'Map data &copy;  OpenStreetMap  contributors, ' +
    ' CC-BY-SA , ' +
    'Imagery ©  Mapbox ',
            key: 'BC9A493B41014CAABB98F0471D759707'
        })
};

var map = L.map('map', {
    center: new L.LatLng(44.9323281,26.0306833),
    zoom: 12.25,
    layers: [tileLayer['Harta'], groupA, groupB, groupC] //change this to determine which ones start loaded on screen
});

//Control on the Top Left that handles the switching between A and B
var overlayMaps = {
    "Scoli": groupA,
    "Muzee": groupB,
    "Parcuri": groupC
};
L.control.layers(tileLayer, overlayMaps, {position:'topleft'}).addTo(map);


texts: Reading and setting label settings in PyQGIS 3?

texts:
How do I read and set label style settings in QGIS 3.0.

For instance, in QGIS 2.18 i had done:

props = {
    "labeling": "pal",
    "labeling/enabled": "true",
    "labeling/fontFamily": "Arial",
    "labeling/fontSize": "12",
    "labeling/fieldName": "Name",
    "labeling/bufferDraw": True,
    "labeling/bufferSize": "1",
    "labeling/bufferColor": "white"
    }
for k, v in props.items():
    my_layer.setCustomProperty(k, v)

How do I do the equivalent in QGIS 3.0?


texts: How to read non-spatial data using PyQgis?

texts:
I'm generating maps with QGis. What I do is generate a bunch of non-spatial data on a set of cities (using VBA), and still using VBA, generate a python code that will automaticaly output general map on the studied territory.

So far, I was importing the data into Qgis with the following code :

lyrData = QgsVectorLayer(u'C:/SomePlace/SomeFile.xlsx', ' ', 'ogr')
QgsMapLayerRegistry.instance().addMapLayer(lyrData)

Which seemed to work, until it started to treat the header of the attributes/columns as data for some territories. Note that this bug is consistent : for some territories it will import header as header and for some it will import header as data...

I tried different things so far : If I import manually (without pyQgis), the bug is still here. I tried renaming the .xlsx file with the name of one that works, but it didn't work.

Another thing I tried is to import the file as an .csv instead of .xslx, but Qgis Cookbook only give exemple to import spatial data csv :

uri = "/some/path/file.csv?delimiter=%s&xField=%s&yField=%s" % (";", "x", "y")
vlayer = QgsVectorLayer(uri, "layer_name_you_like", "delimitedtext")

Or

uri = "file:///some/path/file.csv?delimiter=%s&crs=epsg:4723&wktField=%s" % (";", "shape")

And this, despite the fact that you can manually import CSV files without geometry

Using PyQGIS, is there a way to precise that imported excel files have an header, or is there a way to import csv files without geometry?


texts:
I am somewhat disheartened by my experience with Spatialite databases in QGIS. In theory, it could be what personal geodatabases are for ArcGIS: a consolidated, portable format that keeps projects clean. In reality, I run into problems though, that make it very cumbersome to use Spatialite in QGIS:

Data types do not get assigned properly, most fields get TEXT, even if the CREATE statement stated otherwise. This makes form creation impossible since it depends on the data type of a field.
The GUI accepts Spatialite, however, the output format is almost always a Shapefile, truncating attribute names in the process… Boolean operations between Shapefiles and Spatialite do not work at all…
The functionality to work with other formats is well hidden in the UI. in the Toolbox there are the socalled "QGIS-Geo-Algorithms" that basically do all the Vector functions with other formats as output. However, it is not possible to specify a spatialite table as the target, it just works with spatialite databases.

Is there a good way to use Spatialite in QGIS?


texts: QGIS 2.0 Python API: Layer vector style loaded from file, marker symbology is respected, labels are not

texts:
Im loading a QgsVectorLayer and showing it on map. I load the style i prepared with desktop QGIS 2.0.1 using the layer.loadNamedStyles(path) function in the API and it seems to work, kind of, because the markers are shown on map canvas the way i specified, however the labels are not shown even though the labels styling is defined in the .qml style file (looked at it in notepad). Is there something else i need to do, to display the labels for a vector layer?

Result is the same if i save the style as .sld then load it via layer.loadSldStyle(path)

I think the issue is the rendering engine/settings used. The style files assumes 'pal' rendering in the custom properties, so it could be something with the way Pal related objects are (not) initialized for the layer

This seems to have enable the labels, but it just displays "label" for every feature:

layer.enableLabels(True)

I still think the style i prepared is being ignored for some reason. There must be a simpler/better way to do this!

If i use a qml file, generated with qgis 1.8, labels show up without me doing anything other than loading the layer. Must be 2.0.1 related.

UPDATE: I finally know what the issue is: .qml style file. Despite the fact that it was generated with qgis 2.0.1 it seems that when style file is loaded into my stand alone app, it tries to load the label styling information from the legacy part of the file, which is not filled out:

  <displayfield>name</displayfield>
  <label>0</label>
  <labelattributes>
    <label fieldname="" text="Label"/>
    <family fieldname="" name="MS Shell Dlg 2"/>
    <size fieldname="" units="pt" value="12"/>
    <bold fieldname="" on="0"/>
    <italic fieldname="" on="0"/>
    <underline fieldname="" on="0"/>
    <strikeout fieldname="" on="0"/>
    <color fieldname="" red="0" blue="0" green="0"/>
    <x fieldname=""/>
    <y fieldname=""/>
    <offset x="0" y="0" units="pt" yfieldname="" xfieldname=""/>
    <angle fieldname="" value="0" auto="0"/>
    <alignment fieldname="" value="center"/>
    <buffercolor fieldname="" red="255" blue="255" green="255"/>
    <buffersize fieldname="" units="pt" value="1"/>
    <bufferenabled fieldname="" on=""/>
    <multilineenabled fieldname="" on=""/>
    <selectedonly on=""/>
  </labelattributes>

If i changed the file by hand accordingly, the style now works. But the questions is: how do i force, or properly use custom properties where the style i saved was actually saved?


texts: I'm using QGIS 2.18.3 with the cloud plugin. The map that I've published using the plugin doesn't seem to be retaining the coordinate display setting that I specify in Project-Project Properties-Coordinate Display. Could this be a problem with the plugin, or is there some other place I have to specify the setting to get it to display correctly in a published map?  Other project properties seem to publish just fine.


texts: Proj4J not precise for EPSG:3857 transformations

texts:
I try to convert some coordinates from "EPSG:31469" to "EPSG:3857" and back.

With the heavy framework GeoTools it seems to work pretty precise.

"EPSG:31469" - 5439627.33,          5661628.09
"EPSG:3857"  - 1573655.6648492748,  6636624.730959651
"EPSG:31469" - 5439627.330626475,   5661628.099510074

But I would like to use the lightweight library Proj4j (https://github.com/Proj4J/proj4j) which produces massive deviation:

"EPSG:31469" - 5439627.33,          5661628.09
"EPSG:3857"  - 1573659.9666159092   6603334.624358153
"EPSG:31469" - 5439626.576916553    5661585.491150079

Any ideas why this is the case?

My code:

import org.cts.crs.CRSException;
import org.osgeo.proj4j.CRSFactory;
import org.osgeo.proj4j.CoordinateTransform;
import org.osgeo.proj4j.CoordinateTransformFactory;
import org.osgeo.proj4j.ProjCoordinate;

public class Test
{

    public static void main(String[] args) throws CRSException
    {
        ProjCoordinate coord = new ProjCoordinate(5439627.33, 5661628.09);
        System.out.println(coord);

        ProjCoordinate target = new ProjCoordinate();
        CRSFactory crsFactory = new CRSFactory();
        CoordinateTransformFactory f = new CoordinateTransformFactory();
        CoordinateTransform t;

        t = f.createTransform(crsFactory.createFromName("EPSG:31469"),
            crsFactory.createFromName("EPSG:3857"));

        t.transform(coord, target);

        System.out.println(target);

        ProjCoordinate coord2 = new ProjCoordinate(target.x, target.y);

        t = f.createTransform(crsFactory.createFromName("EPSG:3857"),
            crsFactory.createFromName("EPSG:31469"));

        t.transform(coord2, target);

        System.out.println(target);
        System.out.println(target.x - coord.x);
        System.out.println(target.y - coord.y);

    }

}


texts: QGIS print layout labels size respective to polygon shape

texts:
I am working with a tax lot shapefile with attributes such as ownership, addresses, account numbers, etc. I am trying to create a map with the ownership name of each polygon. My problem is that some of the polygons are huge and some are small, so when I go to my print layout page with those labels turned on, the large polygons have a tiny bit of text in them and the small polygons have text that falls outside of it's boundaries.
I would assume a simple if else statement would fix something like this. Though I am not sure how to do and if else statement for a label nor have any experience using python within QGIS if that is the route for this problem. Could someone help guide me though this process?


texts:
I have different layers with different CRS, I think therefore the layers don't match. But I don't know how to fix this. Here is what I did (QGIS 3.0.1 on Debian Buster).

I connected to OpenStreetMap Tile Server https://a.tile.openstreetmap.org/{z}/{x}/{y}.png to request XYZ Tiles. This is the first layer.

The tile server uses the CRS WGS 84, Pseudo Mercator Projection EPSG 3857.

Now, I tried to built my own map, by doing a quick query with QuickOSM. Key is building, querying Extent of the map canvas (coordinates are 300.000, 6.100.000, Scale is 1:60.000, which is in France, somewhere south of Paris).

This gives me three additional layers, all called OsmQuery, see image.

The additional layers use CRS WGS 84 EPSG:4326. If you look at the map, you see the buildings belonging to Chateau-Landon (south) near La Madeleine-sur-Loing (north). The destinations are about 7km apart.

My thinking is: "The error might result from different projections. However, 7km is a lot. I find it kind of strange that OpenStreetMap uses two different CRS." But I'm quite a newbie...

What can I do to get an exact match of all the layers?

My QGIS version is as follows:


texts:
How do I add a custom coordinate system to FME?

I need to transform my data to EPSG 3765 which is not in Coordinate System Gallery.


texts: Applying css on polylines using Leaflet?

texts:
I am working with the application which uses leaflet api.

Introduction

I needed to draw different types of fences, using decorators i can somewhat apply good visuals to the polylines but not much.

Problem

I was willing to show twisted wires instead of dashes, dots or plain lines and I know the twisted wire line will be an image but can't find help about applying custom css to polylines.

Script Example

         var fence2Icon = L.icon({
                            iconUrl: 'xxxx.png',
                            iconSize: [5, 20]
                            iconAnchor: [5, 18]
                        });

                        // Add coordinate to the polyline
                        var polylineFence2 = new L.Polyline([], { color: 'red' });
                    function fencePlace2(e) {
                        // New marker on coordinate, add it to the map
                new L.Marker(e.latlng, { icon: fence2Icon, draggable: false }).addTo(curr);
                        // Add coordinate to the polyline
               polylineFence2.addLatLng(e.latlng).addTo(curr);
            var decorator = L.polylineDecorator(polylineFence2, {
            patterns:[{offset:5,repeat:'20px',symbol:new L.Symbol.Dash({pixelSize:5})
                         }]
                    }).addTo(curr);
                    }

                L.easyButton('fa-flash', function () {
                            $('.leaflet-container').css('cursor', 'crosshair');
                            map.on('click', fencePlace2);
                            polylineFence2 = new L.Polyline([], { color: 'red' });
                        }).addTo(map);


texts: Data frame to SpatialPolygonsDataFrame with multiple polygons

texts:
I have a data frame with a bunch of polygons, several points per polygon. How can I convert it to a SpatialPolygonsDataFrame object?

Here's my data frame:

> buildings_df
        lon      lat               id
1  96.13597 16.83363 convocation_hall
2  96.13580 16.83315 convocation_hall
3  96.13631 16.83298 convocation_hall
4  96.13656 16.83351 convocation_hall
5  96.13400 16.83104    judson_church
6  96.13392 16.83073    judson_church
7  96.13432 16.83064    judson_church
8  96.13444 16.83085    judson_church
9  96.13596 16.83074      the_library
10 96.13586 16.83044      the_library
11 96.13615 16.83035      the_library
12 96.13626 16.83065      the_library
13 96.13599 16.83073      the_library

This can be reproduced with:

> dput(buildings_df)
structure(list(lon = c(96.135972, 96.1358, 96.136315, 96.136562,
96.133997, 96.133922, 96.134319, 96.134437, 96.135961, 96.135864,
96.136154, 96.136261, 96.135993), lat = c(16.833625, 16.833153,
16.832978, 16.833512, 16.831038, 16.830729, 16.830637, 16.830853,
16.83074, 16.830442, 16.83035, 16.830647, 16.830729), id = structure(c(1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L), .Label = c("convocation_hall",
"judson_church", "the_library"), class = "factor")), .Names = c("lon",
"lat", "id"), row.names = c(NA, -13L), class = "data.frame")

From this data.frame, I then convert it to a list, then to polygons, then to a SpatialPolygons object

# make a list
buildings_list <- split(buildings_df, buildings_df$id)
# only want lon-lats in the list, not the names
buildings_list <- lapply(buildings_list, function(x) { x["id"] <- NULL; x })

#  make data.frame into spatial polygon, cf. http://jwhollister.com/iale_open_science/2015/07/05/03-Spatial-Data-In-R/

# create SpatialPolygons Object, convert coords to polygon
ps <- sapply(buildings_list, Polygon)

# add id variable
p1 <- Polygons(ps, ID = 1)

# create SpatialPolygons object
my_spatial_polys <- SpatialPolygons(list(p1), proj4string = CRS("+proj=longlat +datum=WGS84") )

# let's see them
plot(my_spatial_polys)

So far so good:

And each of the polygons seems to be in there...

> str(my_spatial_polys)
Formal class 'SpatialPolygons' [package "sp"] with 4 slots
  ..@ polygons   :List of 1
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 3
  .. .. .. .. ..$ convocation_hall:Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 16.8
  .. .. .. .. .. .. ..@ area   : num 3.08e-07
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. .. ..$ judson_church   :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 43.3
  .. .. .. .. .. .. ..@ area   : num 0.0204
  .. .. .. .. .. .. ..@ hole   : logi TRUE
  .. .. .. .. .. .. ..@ ringDir: int -1
  .. .. .. .. .. .. ..@ coords : num [1:5, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. .. ..$ the_library     :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 96.1 43.3
  .. .. .. .. .. .. ..@ area   : num 0.00512
  .. .. .. .. .. .. ..@ hole   : logi TRUE
  .. .. .. .. .. .. ..@ ringDir: int -1
  .. .. .. .. .. .. ..@ coords : num [1:6, 1:2] 96.1 96.1 96.1 96.1 96.1 ...
  .. .. .. ..@ plotOrder: int [1:3] 1 2 3
  .. .. .. ..@ labpt    : num [1:2] 96.1 16.8
  .. .. .. ..@ ID       : chr "1"
  .. .. .. ..@ area     : num 0.0302
  ..@ plotOrder  : int 1
  ..@ bbox       : num [1:2, 1:2] 96.1 16.8 96.1 16.8
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0"

But they don't seem to be grouped in the way I'm expecting. When I try to set three ID variables, I get an error. How can I get the SpatialPolygons to have three polygons, each with its own ID?


texts: Show a single value related to a whole layer in QGIS

texts: I'm developing a maximum flow plugin for QGIS.  I would like to show a calculation result (i.e., the maximum flow) on the map.  The result is not a property of a single feature.  It is a property of a whole layer.  What is the best way to accomplish this?