texts:
After I have created a new QGIS project, the project CRS is set automatically to WGS84 (EPSG:4326 - WGS 84 - Geographic). The layer (.shp) which I loaded to the project has already a predefined CRS (also WGS84). I want to change the CRS to a projected CRS to make some calculations so I set the appropriate CRS both for the layer and for the project (EPSG:32632 - WGS 84 / UTM zone 32N – Projected). In the project properties (General Tab) WGS84 is set as the ellipsoid for distance and area calculations.

Although both the project and the layer(s) have exactly the same CRS, the layer is displayed in a geographical CRS (WGS84):

Changing the layer CRS back to WGS84 I actually get the correct display (as a projected CRS - UTM32N) although the layer CRS is wrong (WGS84):

To check the correct position I also loaded OSM as basemap.
Using WGS84 I get the following result (exptected for UTM...):

project: WGS84

layer: WGS84

OSM basemap: WGS84

ellipsoid: WGS84

while by:

project: UTM32N

layer: UTM32N

OSM basemap: UTM32N

ellipsoid: WGS84 or GRS80

I am landing in the western coast of Middle Africa.

This is really annoying. Do I misunderstand anything generally regarding CRS settings in QGIS or what is the problem?


texts: Using Image exportToDrive in Google Earth Engine?

texts:
I want to export image to drive for my roi. For this i have specified a geometry(roi). But exported map doesn't cover the area specified and it's not showing the color palette. Is there any possible way to export a colored image (which shows in my layer) to Drive? Why it's not covering the region I have specified.
Problem is in export.To.Drive section. But all the codes are given to check error that i did.

Code was:

var region = ee.FeatureCollection('users/karmakersourav/shape')
  .filter(ee.Filter.eq('NAME_1', 'Chittagong'));
var addNDVI = function(image) {
  var ndvi= image.normalizedDifference(['B4', 'B3']).rename('NDVI');
  return image.addBands(ndvi);
};

var cloudMask = function(image) {
  var clouds = ee.Algorithms.Landsat.simpleCloudScore(image).select(['cloud']);
  return image.updateMask(clouds.lt(25));
};
// Select the land mask.
var landMask = function(image) {
var datamask = hansenImage.select('datamask');
// Create a binary mask.
var mask = datamask.eq(1);
return image.updateMask(mask);
};
// Load a Landsat collection, apply the NDVI and cloud/water masking functions over it.
var collection1 = ee.ImageCollection('LANDSAT/LE07/C01/T1_RT_TOA')
  .filterBounds(region)
  .filterDate('2008-01-01','2008-12-31')
  .map(addNDVI)
  .map(cloudMask)
  .map(landMask)
  .select('NDVI') ;
print(collection1);

//reducer
var annualNDVI = collection1.reduce(ee.Reducer.median());
var region_NDVI = annualNDVI.reduceRegions({
  collection: region,
  reducer: ee.Reducer.median(),
  scale: 30 // the resolution of the dataset
});

var ndviParams = {min: -1, max: 1, palette: ['blue', 'white','green']};
var boundParams = {palette: ['black']};
Map.addLayer(annualNDVI.clip(region), ndviParams, 'NDVI', 1);
Map.addLayer(region, boundParams, 'LA', 1);

//here's the problem

Export.image.toDrive({
 image: annualNDVI,
 description: 'CTG_NDVi',
 scale: 30,
 region: roi,
 fileFormat: 'GeoTIFF',
 maxPixels: 1e12,
 formatOptions: {
    cloudOptimized: true
  }
});

My imports was:

var image = ee.Image("UMD/hansen/global_forest_change_2015_v1_3"),
    region = ee.FeatureCollection("users/karmakersourav/shape"),
    hansenImage = ee.Image("UMD/hansen/global_forest_change_2015"),
    roi = /* color: #98ff00 */ee.Geometry.Polygon(
        [[[91.26571433292804, 22.382447200369512],
          [91.70516745792804, 22.24015355729955],
          [91.81503073917804, 21.50608889294255],
          [91.99081198917804, 21.158146277513005],
          [92.36434714542804, 20.74775339022756],
          [92.38631980167804, 21.362918378188724],
          [92.68868864595584, 21.285515736407884],
          [92.66671598970584, 22.12251851402064],
          [92.42501677095584, 22.964696114885136],
          [92.43600309908084, 23.237540932485988],
          [92.31515348970584, 23.620600408965007],
          [92.21627653658084, 23.781555321668222],
          [92.04049528658084, 23.691042604440895],
          [91.90865934908084, 23.731278240693108],
          [91.89767302095584, 23.47960242153599],
          [91.75485075533084, 23.28800680321665],
          [91.77682341158084, 23.08602889002461],
          [91.61202848970584, 22.99503954262493],
          [91.46920622408084, 23.277915158047364],
          [91.34835661470584, 23.308187798303887],
          [91.22750700533084, 23.600467103459394],
          [91.28243864595584, 23.872005072934794],
          [91.32638395845584, 24.082808607079368],
          [91.24947966158084, 24.263221977268092],
          [91.06271208345584, 24.223152164577193],
          [90.73312223970584, 23.88205115062075],
          [90.63424528658084, 23.660857737733508],
          [90.54635466158084, 23.328365731356676],
          [90.65621794283084, 22.83312931280329],
          [90.98580778658084, 22.437667249990497],
          [90.98580778658084, 22.071621721722586]]]);


texts: Is there a better way to convert Raster to Point data? Don't need one point for every cell

texts:
I need to convert raster data to point data so I can use the points for interpolation via kriging. My problem is that my resolution is 5 ft and I'm working on a whole city.  The data I am converting to points is grouped by zone with No Data values between; most/all of the points in each zone have the same value. I need the high resolution because some zones are linear, just 2-3 cells wide and about 20 cells long, while others are huge behemoths.  I need to capture the smaller zones but I don't need 2000 points in the larger zones; because the values are all the same, I might need 10 or 20.  I experimented with lowering the resolution to 50 and 100 ft but I lose much of the smaller zones. My raster is floating data so converting to polygons to do any other vector stuff isn't possible, I don't think.

Am I doomed to have a billion points and manually remove ones from the big data chunks?

What kind of processing time am I looking at if I do 5 ft Raster to Point and have one point per cell and I interpolate those?  It already took several hours to interpolate between 10 points via kriging over the whole area.

I would not like to sacrifice my data resolution for the small zones if I can help it.  They are very important and represent the variability I need to capture.


texts:
I'm using openlayers 4 to render rivers as shown on the screenshot below (map zoom level 7) and would like to display river names along the lines. The river features are openstreetmap waterway=river ways served by geoserver to a vectore layer.
The features are styled like this, but the lines are so curvy that the names are almost never rendered.
var styles = {
'river_lines': new ol.style.Style({
    stroke: new ol.style.Stroke({
        color: 'blue',
        width: 1.25
    })
}),
'river_labels': new ol.style.Style({
    text: new ol.style.Text({
        text: '',
        font: '15px sans-serif',
        fill: new ol.style.Fill({
            color: 'white'
        }),
        placement: 'line',
        maxAngle: 0.15
    })
})
};

I've tried replace the geometry of the label style with a simplified one but it did not change anything.
styles.river_labels.setGeometry(feature.getGeometry().simplify(0.9));

Any ideas how to solve this?


texts: Using PostGIS functionalities in Microsoft SQL Server?

texts:
I'm using Postgresql 8.4, Postgis 1.5 for My GIS works.
Planning to change my Database to Microsoft SQL Server 2012.

Can we use all the Postgis functionalities in Microsoft SQL Server 2012?


texts: How to create updateable/inseartable view on QGIS with a sequence id?

texts:
I have a huge point dataset that takes a long time to insert new stuff when using QGIS

This table has an id column of type serial

To go around this problem I've created an empty view that has the same exact fields as its original table. at first tried limit 0 but it didn't let it update so I limited the table with a where condition that never gets fulfilled, so i end up with an empty table

Problem now is that when inserting new features, they don't get the id of the original table, and new inserts on the original table need the id to activate other triggers

How to make it use the same sequence as the original table and enable people to insert new features on it?


texts: How to properly label colliding labels in QGIS 3.x

texts:
Is there a way to properly show the labels when "Show all labels for this layer (including colliding labels)" is selected from label properties? I was looking for adding callout lines for the labels but I think it requires the labels to be manually displaced. As there's more than 1000 points it will be too time consuming to do this. Is there an easier way to do this?

It should be sufficient if the points are clustered and each point in that cluster can be labelled as shown below.


texts: Creating DEM from point data in CSV files using QGIS?

texts: how can we create a DEM in QGIS from point data in CSV files? A similar problem was discussed here, using a shp-file. Our CSV's are too big, so that dxf2shp crashes. So we are looking for another way to do this.


texts: Is there a bug in the Build Virtual Raster tool (QGIS 3.4.7) when using nodata value = -9999?

texts:
When I configure the Nodata Value argument in the Build Virtual Raster tool (QGIS 3.4.7) equal to -9999 I have obtained an unexpected result, as shown in the figure below:

On the other hand, if I adopt Nodata Value = 0 the result is better, but without the representation of pixels with values = 0:

Raster layer cuts were done with the GDAL Clip Raster by Mask Layer tool (Nodata value = -9999, Match the extent of the mask layer, Profile: default, Output data type: Float64).

Is this a bug in the QGIS 3.4.7 Build Virtual Raster tool?


texts: Debugging QGIS plugins using Eclipse IDE in Ubuntu?

texts: I am using QGIS 2.14.1 in Ubuntu 15.04. I am developing QGIS plugin using Eclipse 3.8 IDE installed. I would like to configure pyqgis and debug the code using eclipse. Is there any documentation for configuring settings for developing plugins using eclipse for Ubuntu 15.04?


texts: Where to download models to run in ModelBuilder?

texts:
I am wondering if there is a website where I can download models to run in ModelBuilder?

A place where people post and share different geoprocessing models to be run inside of ModelBuilder.


texts: Getting band names using ArcPy in ArcGIS Desktop 10.x

texts:
How do I get a list of band names using ArcPy in ArcGIS Desktop 10.x?
In ArcGIS Pro, I can use Raster('myRaster').bandNames, but this is not an option in ArcGIS 10.x and I cannot seem to figure out how to access band names.


texts: Releasing ArcMap memory by removing *.sr.lock files from shapefiles

texts:
Is there a method of removing *.sr.lock files from shapefiles to release memory consumption without closing out of an ArcMap application?  I am attempting to run code that will cycle through 200+ folders filled with shapefiles and perform various actions on them.  I close the current map document and open a new one after each folder.  However, this is causing ArcMap to build up memory and eventually crash after around 100 folders are processed.  A big part of this problem (I would suspect) is the existence of lock files on each shapefile that exists. I have tried using a ComReleaser and Marshal.ReleaseComObject, but nothing seems to get rid of these lock files aside from completely closing out of ArcMap, which is not an option.  Is this just the way ArcMap handles shapefiles?
The following is a small portion of my code that creates the lock files:
// Loop through the selected shapefiles and load them into the map

pGxObjects.Reset();

while ((pGxDataset = (IGxDataset)pGxObjects.Next()) != null)
{

    pProDlg.Description = "Loading selected shapefile(s)...";
    pStepPro.Message = pGxDataset.DatasetName.Name;

    IFeatureLayer pFeatLayer = new FeatureLayer();
    pFeatLayer.FeatureClass = (IFeatureClass)pGxDataset.Dataset;
    pFeatLayer.Name = pFeatLayer.FeatureClass.AliasName;

    pMxDoc.FocusMap.AddLayer(pFeatLayer);

    // Check if the cancel button was pressed. If so, then stop the process.
    boolCont = pTrkCan.Continue();
    if (!boolCont)
    {
        pProDlg.HideDialog();
        return null;
    }

    pStepPro.Step();
}


texts: Creating legend in separate page using ArcMap?

texts:
Is it possible in ArcMap 10.2.2. to create a layout that shows the legend of the map on a separate page? A.k.a have two pages in layout view?

So first page: full view of the map, second page: legend.

Later on I would like to create an ArcGIS Server printing service from this layout, and use it in WebApp Builder for ArcGIS, so that's why it is important to integrate both of it in one layout.


texts: Detecting overlap for buffered lines using ArcGIS Desktop?

texts:
We have almost 50000+ lines which we want to remove a line if it is too close to the other.

So we tried to do a buffer for two line, then check if the buffer area are overlapped.

I wonder if there is a efficient idea to make it?

Since at my first thought, we need almost 50000*(50000-1) times detect for two buffered area, it will take a long long time plus the buffer time for a line. :(

We use ArcGIS, but we prefer to make it by program rather than through the UI.

Also, we will give a certain distance of meters or pixels to tell the program if the two line are too close to each other.


texts:
I am reopening an issue about:

  RuntimeError: ERROR 000824: The tool is not licensed

Since I could not find the answer to my problem in previous issues.

I am running a script.py in Anaconda. The scripts requires arcpy that I can import to python.

When I run the script, everything seems alright until the end where I get the following error:

File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4483, in Times
    in_raster_or_constant2)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Utils.py", line 53, in swapper
    result = wrapper(*args, **kwargs)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.5\ArcPy\arcpy\sa\Functions.py", line 4480, in Wrapper
    ["Times", in_raster_or_constant1, in_raster_or_constant2])
RuntimeError: ERROR 000824: The tool is not licensed.

I have checked my extensions in ArcMap and all of them are ticked. In addition, I have a license for ArcGIS that I had previously authorized.

I cannot figure it out and I don't know which tools are not licensed.


texts: Why are services directory and arcgis manager not accessible?

texts:
I installed ArcGIS Server 10.0 for Java with its post install and its SP 5.

When I try to access the services directory, I get a 404 error from the web server. I don't get that error when I try to get to the arcgis manager web interface.

After displaying the arcgis manager web interface, when I try to login, I get an error: Login information is incorrect. The login used is arcgismanager that is member of agsadmin. The password used is the one specified in the post install.

Do you have any hint?


texts: How to export OSM maps given the bounds?

texts:
I have a map application where I use all kinds of OSM maps as background.

http://{a-c}.tile.openstreetmap.org/{z}/{x}/{y}.png
http://{a-c}.tile.opencyclemap.org/cycle/{z}/{x}/{y}.png
http://{a-c}.basemaps.cartocdn.com/light_all/{z}/{x}/{y}.png
http://{a-c}.tile.stamen.com/watercolor/{z}/{x}/{y}.png

Is there a way to extract a part of the OSM map given the bounds?

I looked here and I tried this:

http://a.basemaps.cartocdn.com/api/0.6/map?bbox=-10,-10,10.10

But it didn't work.

I want to extract the part of the map that is within those bounds and get it returned as an image.


texts: ArcGIS 10.2 Custom Geocoder outputs results with no geometry

texts:
I have a custom general single field geocoder that reads a meter number field attached to premise points. This meter number field is the only link between our premise points and corrosion data tables with location information.

I created this geocoder so that our corrosion team can geo-locate keypoints based on meter numbers. The geocoder works fine (the correct coordinates are generated and the correct meter numbers are matched) however, the geocoder does not output any geometry in the results, just data. I have verified the projected coordinate systems and the source data of the geocoder.

Could something be wrong with the parameters of my geocoder? Our temporary solution to this is to extract the geocode result table and simply display XY data.


texts: What is the reason to use GeoServer to provide access to vector data?

texts:
I am building a website that will present spatial data (Leaflet) and I plan to include functionalities like: add, move, delete features (i.e. by clicking on the map or on the feature) I would like to use PostGIS as DB.
Should I use a GeoServer and if so why?
I was planning to store vector data and request them as GeoJSON from DB. The approximate amount of records (10 000 - 100 000).


texts: Using File GDB API .NET Wrapper in ASP.NET Web Service

texts:
I am trying to use the File GDB API (1.2) .NET wrapper in an asp.net web service with no luck.  The website will not load spitting out error:  "Could not load file or assembly 'Esri.FileGDBAPI.DLL' or one of its dependencies. The specified module could not be found."

The .NET Wrapper works perfectly fine in a windows form application from my machine.  Just to rule out the obvious:

Esri.GileGDBAPI.dll, FileGDBAPI.dll, and FileGDBAPID.dll all are sitting in the bin of the website (along with OpenTK dll's)
I have the prerequisite OpenTK installed
Latest VC Runtime distributable installed.
Website is targeting .NET 4.0

The fact that it runs fine in a windows form application gives me confidence that I have everything I need on my machine.

Does anyone have any other ideas and/or has anyone successfully used the .net wrapper in a .net webservice?

UPDATE:  With the help of fuslogvw in the post I marked as the answer, I was able to to get the web service working after the directory containing the FileGDBAPI.dll and FileGDBAPID.dll was added to the 'Path' variable.


texts:
I'm trying to connect data from two different layers to combine them in one new layer and continue working with them.
Let me explain that a bit more detailed: The first layer contains data (and borders) from counties whereas the other layer contains operating areas of businesses units (with different borders than the counties).
My aim is to compare revenues from the operating areas with county data (such as purchasing power). Is there such a function?

I'm new to QGIS.

I added a sketch to illustrate the problem. The operating areas are sometimes within two counties.


texts: Seeking maps to use on iPhone/Android app

texts:
How can I get data of the earth map on the iOS sdk "Swift" containing data altitude too?
I would like to use all world map, gps position (but I think that's not a matter) and altitude of all places.
Is there an option that has less computational cost?
I've tried only Google Maps but I've read that if I will public my app release and only a few people will use it I will have to pay because it will use a lot of data.


texts: What is the reason for the very poor performance of ESRI's file geodatabases on network drives?

texts:
I am asking this question because I am at a loss and dont have a good enough understanding of all the underlying operations going on 'under the hood' of Windows/Arcgis.

What is the root cause of fgdb performance issues on network drives?

I have done some testing regarding performance of fgdb's vs shapefiles, network vs local drive.

B: drive is my local computer drive.
K: drive is the local network drive for geodata.

Method: I use 5 point feature classes with random generated points, all within the same extent. The feature classes have respectively 100 000, 10 000, 1000, 100 and 10 points. The left side of the diagram shows time used while transferring higher number of points, while the right side shows time used while exporting feature classes with 10 points (all duplicates of the 10 point file). The Y-axis is time in seconds.

Note that all results on my local computer drive are pretty much flat compared to transfers on the network.

Monitoring the task manager on the network tab on my computer, i can see that there are large periods of network 'silence' followed by a short burst of activity (1 featureclass), then followed by anoth 40'ish seconds of nothing.I am asuming it is because of schema locks, But I would like to know what causes it and what can be done to prevent it.

Backround information:
I work within a very large government organization (~1300'ish ArgGis licenses) with several regional offices and a small'ish central IT group. We are still working in ArcGIS 9.3.1, but are going to upgrade to 10.3 in the fall. We are not a GIS-dedicated organisation, GIS is only a tool for supporting the decisionmakers here. Only my colleague and I are GIS-engineers in regional office (500 employees, 70 Arcgis licenses).
We have a central database (SDE on Windows Server 2008 R2 that contains most data on a national level) and a local (physical server) that is dedicated for geodata within each regional office. Certain datasets (large raster datasets, DEMs) are divided on a national level and also stored on the local server in the regional office.
The infrastructure is set up by IT. Geodata is only a small part of this, so the geodata considerations do not always take priority with the entirety of the IT system (i.e. turning off opportunistic locking is not an option).
Due to our large organization and legal requirements, working on the local computer drive should be avoided due to lack of backup. In theory, all geodatawork should be done on the geodata network drive (local in-house server).
When we upgrade to ArcGIS 10.3, with the easy option of 'let this be my home directory for this MXD' option, I would very much like to use the filegeodatabases to help keep all data organised, as there are several new MXD per week with plenty of scratch data floating around everywhere. However, with the performance issues that I have found, It feels we will still be stuck with shapefiles and a hard-to-keep-organized data environment.


texts: python GDAL 'GetDriverByName', argument 1 of type 'char const *'

texts:
i want to use python GDAL to convert some rasters files but i take error in the create new tiff image and specific in the line driver= gdal.GetDriverByName('GTiff')

update :

python version 2.7 64 bit

GDAL version 2020100

os: windows

i use it and before this version and all work fine
code :

inputfile = 'input.tif'
input = gdal.Open(inputfile )
inputProj = input.GetProjection()
inputTrans = input.GetGeoTransform()

referencefile = 'ref.tif'
reference = gdal.Open(referencefile )
referenceProj = reference.GetProjection()
referenceTrans = reference.GetGeoTransform()
bandreference = reference.GetRasterBand(1)
x = reference.RasterXSize
y = reference.RasterYSize

outputfile = tempfile.NamedTemporaryFile(suffix='.tif').name
driver= gdal.GetDriverByName('GTiff')
output = driver.Create(outputfile,x,y,1,bandreference.DataType)
output.SetGeoTransform(referenceTrans)
output.SetProjection(referenceProj)

gdal.ReprojectImage(input,output,inputProj,referenceProj,gdalconst.GRA_Bilinear)

error :

    return _gdal.GetDriverByName(*args)
TypeError: in method 'GetDriverByName', argument 1 of type 'char const *'

update 2:

that error show when i use this limes:

from __future__ import absolute_import
from __future__ import unicode_literals

but that i need it because i use python 2.7

any idea how to fix it ?


texts: Add-Ins Category missing from Customize / Commands / Categories - Can't Use Add-Ins?

texts:
I've created a Add-In button that says "Hello World" when clicked.
I want to add button to a toolbar.

I've tried creating my own "Category" in the Add-In wizard, no success.
I've also accepted the default "Category" called "Add-Ins", no success.

The Add-In appears to load but is not assigned to a visible "Category".
So, I can't use it.

I've run the Add-In via Visual Studio, Installed via the ESRI ArcGIS Add-In Installation Utility, and loaded the file via the Customize / Customize Mode / Commands.

What is the issue?
Do I need to enable Add-Ins or turn on an extension?


texts: Convert list of coordinates to OGRGeometry or WKT


texts: Copying feature classes from gdb to different folder locations in Python?

texts:
I am very new to Python, and I'm trying to create a function where I'm copying over feature classes from a gdb to different folder locations as shapefiles. The script needs to be reading from a textfile that lists the different path locations. Here is a screenshot of the textfile:

So if the feature class from the gdb is "apd" it needs to copy over to "\DataSets\APD\apd.shp", if it's "comps", it needs to be copied over to "\DataSets\Completion\comps.shp", and so on. The function needs to have two parameters (the input and the output). The input being the feature class and the output being the destination path.

How do I define my destination parameter so that a particular feature class needs to be copied over to that specific folder?

Sadly, my code is looking something like this:

import arcpy
import os

arcpy.env.workspace = 'GOM3.gdb\\NAD1927'
featureclasses = arcpy.ListFeatureClasses()

def CopyFeature (fc, dest):
    for fc in featureclasses:
        if arcpy.Exists(fc):
                fname="X:\GOM-Cubed\DataSets\ShapefileList.txt"
                with open(fname) as f:
                    content = f.readlines()
                    for dest in content:
                        arcpy.CopyFeatures_management(fc, dest)
                    fname.close()

Eventually I want to add this script to ModelBuilder. That way the user can define the input parameter.


texts:
I have coordinates in WGS 84 34N and I am trying to add the points to ArcGIS 10.2. This should be no problem as I have done it many times but I cannot figure our why it is not working now.

In the Add XY data window I choose the coordinate system WGS 84 34N and an event layer is created. I cannot see the points even if I select a row in the table and click "Zoom to". It zooms to a certain extent which is not correct. There is no point visible in the middle.

What am I doing wrong? I have tried to swap X and Y but it doesn't help. I used an online converter in order to check manually that my coordinates really are located there where they should be, and it is fine.


texts: Get pixel value from a raster in a xy table

texts:
I'm trying to find a way to convert a raster (in my case it's a Landsat satellite image) into a table with 3 columns which would represent: 1) x coordinates (the so called "xllcorners" in ASCII terms); 2) y coordinates (the so called "yllcorners"); 3) the relative pixel values.
As I'm working with 8-bit pixel-depth raster, my table should look like this at the end:

   X       Y      Pixel
493185  4989885     0
493185  4989915    127
   .       .        .
   .       .        .
   .       .        .

I tried to convert the raster into ASCII format and it gave me a grid with 7911 columns and 7211 rows. No problem so far, but since I want to organize the entire grid as it fits my 3 columns requirement, I haven't been able to import it to excel (where I'd like to do my calculations) properely (namely, it seems replacing and cutting off the spaces, inserting commas and so on didn't work so fine as there's always something different than in the "txt" ASCII grid, and it's not a squared grid neither, which complicates the things).
I'm using both ArcGIS 10.1 and QGIS 1.8 if it could be of interest.


texts:
I'm new to OGR module and I can't figure out how to save new layer to my HDD. I allready have shapefile with lines and buffer and I want to clip lines with buffer.

from osgeo import ogr
os.chdir(r"workspace")

# Load shp with lines
driver_lines = ogr.GetDriverByName("ESRI Shapefile")
datasource_lines = driver_lines .Open("lines.shp", 1)
layer_lines = datasource_lines.GetLayer()

# Load shp with buffer
driver_buffer = ogr.GetDriverByName("ESRI Shapefile")
datasource_buffer = driver_buffer .Open("buffer.shp", 1)
layer_buffer = datasource_buffer.GetLayer()

# Create output shapefile
driver_clip = ogr.GetDriverByName("ESRI Shapefile")
datasource_clip = driver_clip.CreateDataSource("clip.shp")
layer_clip = datasource_clip.CreateLayer("clip", geom_type = ogr.wkbLineString)

# Clipping
layer_lines.Clip(layer_buffer, layer_clip)

So I have a layer with clipped lines, but I don't know how to physicaly save it to disc.


texts: Linking Excel spreadsheet property records to cadastral basemap in QGIS using ERF numbers as the unique identifier

texts: How does one link land parcel polygons in QGIS to an Excel spreadsheet of property records, using ERF numbers as the unique identifiers?


texts: Draw Bridges at intersections with correct angle with PyQGIS?

texts:
I have a road layer and a river layer. I want to draw bridges at the intersections with correct orientation of the bridge (parallel to the road).

I have found out that you can use the "Vector->Analysis Tools->Line Intersections" to identify the intersections. But how do I identify the angle of the road layer at this exact point and automatically roate it acordingly?

I have found this rotate point along line layer, they suggested to use the following code:

from math import atan2

# get layers
bridges = QgsMapLayerRegistry.instance().mapLayersByName('bridges')[0]
trails = QgsMapLayerRegistry.instance().mapLayersByName('trails')[0]

# get fieldindex
fni = bridges.fieldNameIndex('angle')

# initializations
tolerance = 0.01  # search tolerance
bridges.startEditing()

# loop over all bridges
for bridge in bridges.getFeatures():
    x = bridge.geometry().asPoint().x()
    y = bridge.geometry().asPoint().y()

    # get the rectangular search area
    searchRect = QgsRectangle(x - tolerance, y - tolerance,  x + tolerance, y + tolerance)

    # find trails
    for trail in trails.getFeatures(QgsFeatureRequest().setFilterRect(searchRect)):
        # get the nearest vertex on trail and the one before and after
        pnt, v, b, a, d = trail.geometry().closestVertex(bridge.geometry().asPoint())
        p1 = trail.geometry().vertexAt(v)
        # when vertex before exists look back, otherwise look forward
        if v>-1 and b>-1:
            p2 = trail.geometry().vertexAt(b)
        elif v>-1 and a>-1:
            p2 = trail.geometry().vertexAt(a)

        # calculate azimuth
        angle = atan2(p2.x() - p1.x(), p2.y() - p1.y()) / 0.017453
        bridge[fni] = angle
        bridges.updateFeature(bridge)

# save changes and stop editing
bridges.commitChanges()

I followed the steps however I get an error for line 15 like so

File "<input>", line 15
bridges.commitChanges()
      ^
SyntaxError: invalid syntax

Line 15 is:

pnt, v, b, a, d = trail.geometry().closestVertex(bridge.geometry().asPoint())


texts: Merging multiple overlapping lines using Python?

texts:
I have a shapefile with different line features, some of which overlap. Is there a way to easily merge all overlapping lines into one, new feature using Python?

I've seen some solutions for overlapping polygons, but have had a lot more trouble finding examples using lines.


texts: Why value relation doesn't work when I delete and connect reference layer?

texts:
I'm trying to connect reference table (domain) with data(text) to my vector layer (id's) using value relation. I have database in Oracle connected. It works fine, but when I delete the reference layer and connect it again, everything is deleted.

Do you know how to make the connection permanent? Or am I doing anything wrong when I am connecting to the database? Or there is any other function to do this connection?


texts: Geoprocessing Widget in ArcGIS online Error: Using symbology layer

texts:
I have created a model in modelbuilder that takes an interactive feature set as a parameter so that the user of the model can draw a polygon "on the fly". I take the interactive feature set and use the Copy Features tool so that I can make the "in_memory" feature set a new feature layer. I have 3 different symbology layers (red, green, yellow) located on the local hard drive of my computer. I use the ApplySymbologyfromLayer_management tool to set the symbology based on a simple if else statement.

This process works with zero problems in modelbuilder in ArcMap. However, when i publish it as a gp service and use it in the geoprocessing widget of ArcGIS Online Web Appbuilder. The symbology layers i have used cannot be found.

Now I know the reason is because the symbology layer is located on my computer but my question is...

Is there a way to use a symbology layer in ArcGIS online web appbuilder for geoprocessing? If that is not possible, what are the alternatives to changing the feature layer to a specific color (i.e red, green, yellow) using arcpy?

I have seen cartographic represenation as a possible workaround however the documentation is kinda confusing and does not really show me how the Rule IDs change the symbology.  I am using ArcMap 10.5 and i do not have Web appbuilder (dev edition).

Here is the error I got from the gp widget in ArcGIS online Web Appbuilder:
ERROR 000732: Symbology Layer: Dataset C:\Documents\SymbologyLayer\yellow.lyr does not exist or is not supported Failed to execute (ApplySymbologyFromLayer).

Here is the code:

def getSymbology(userArea, category):

arcpy.env.overwriteOutput = True
#Copy in memory feature set to in memory copied features
copiedArea =r'in_memory\copiedArea'
arcpy.CopyFeatures_management(userArea, copiedArea)

#Set name of new users feature layer based on category type
if category == "Open Category":
    userLayer = "Open"
elif category == "Specific Category":
   userLayer = "Specific"
else:
  userLayer = 'Certified'

arcpy.MakeFeatureLayer_management(copiedArea, userLayer)

# Set symbology layers
openLayer = r"C:\Documents\SymbologyLayer\green.lyr"

specificLayer = r"C:\Documents\SymbologyLayer\yellow.lyr"

certifiedLayer = r"C:\Documents\SymbologyLayer\red.lyr"

#If category matches the string then apply symbology
if category == "Open Category":
    arcpy.ApplySymbologyFromLayer_management(userLayer, openLayer)
elif category == "Specific Category":
    arcpy.ApplySymbologyFromLayer_management(userLayer, specificLayer)
else:
    arcpy.ApplySymbologyFromLayer_management(userLayer, certifiedLayer)


texts:
i have a problem with some geometry.
When i run the check geometry validity i get an error coming back that the rings are not closed, so recangular forms that i have surveyed are showing up as triangles. When I put the layer in editing mode, the rectangles show their rectangular forms but go back to triangles when i turn off the editing, not sure if this is caused by the codelist on my total station.
a solution i found was to use polygons to lines and then lines to polygons and replace my original polygons with these new ones.
This can be quite time consuming so was wondering if anyone has a quicker fix for this or if there is a way of preventing this.
The problem only occurs with polygons (i use a codelist created on a Leica TC1200)


texts: Problem with MaxScaleDenominator in GeoServer

texts:
I'm using GeoServer version 2.7.6, and I had to publish a WMS service, I tried to establish a limit in the scale (with <se:MaxScaleDenominator>5001.0</se:MaxScaleDenominator>) but when I view the layer, it lets me see it on all scales without limit. That´s to say that MaxScaleDenominator is not working for me.

Please I need help.

<?xml version="1.0" encoding="UTF-8"?>
<StyledLayerDescriptor xmlns="http://www.opengis.net/sld"
   xmlns:ogc="http://www.opengis.net/ogc"
   xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance"
   version="1.1.0" xmlns:xlink="http://www.w3.org/1999/xlink"
   xsi:schemaLocation="http://www.opengis.net/sld   http://schemas.opengis.net/sld/1.1.0/StyledLayerDescriptor.xsd"
   xmlns:se="http://www.opengis.net/se">
  <NamedLayer>
    <se:Name>Cartografia</se:Name>
    <UserStyle>
      <se:Name>Cartografia</se:Name>
      <se:FeatureTypeStyle>
        <se:Rule>
          <se:Name>Single symbol</se:Name>
          <se:MaxScaleDenominator>5001.0</se:MaxScaleDenominator>
          <se:LineSymbolizer>
            <se:Stroke>
              <se:SvgParameter name="stroke">#8c8f8c</se:SvgParameter>
              <se:SvgParameter name="stroke-width">0.1</se:SvgParameter>
              <se:SvgParameter name="stroke-linejoin">bevel</se:SvgParameter>
              <se:SvgParameter name="stroke-linecap">square</se:SvgParameter>
            </se:Stroke>
          </se:LineSymbolizer>
        </se:Rule>
      </se:FeatureTypeStyle>
    </UserStyle>
  </NamedLayer>
</StyledLayerDescriptor>


texts: Cut/Clip/Split a polygon with another polygon from same layer

texts:
To not have overlapping polygons/features, I normally  use one or more polygons to clip/shape other polygons. This is often much faster than tracing when adding polygons etc.
I cannot find this functionality in QGIS.

The highlighted polygon is overlapped by the adjacent (west/left) polygon with the straight edge. I want to cut the adjacent polygon with the highlighted.


texts:
When running my buffer over my points, I am losing all columns into my table.
Here is the query:

SELECT st_buffer(geom, 100)
INTO gis_osm_pois_buf
FROM gis_osm_pois;

Is there any way to keep my attributes in the new table?


texts: 'Scale' value in slope percent map in QGIS?

texts:
What value should be used for scale while calculating slope percent in QGIS?

Some documents suggest is as 111120 and some as 365130.

I am using ESPG 4326 projection and do I need to convert this value to some other value using lat like given in this link:

How to properly get the slope from an ASTER DEM in QGIS?


texts: Calculate shortest distance from point to polygon according angle table

texts: I have a point in polygon. I want to have new polyline shapefile, which will show the shortest distance from the point to polygon edge, according to given angles in the table.


texts: What is the definition for a geographic coordinate system with east/west bounds of 0 to 360

texts: I need the coordinate system reference or a projection definition that uses geographic (i.e. lat/lon) coordinates but instead of east/west bounds being defined as 180 to -180, I need them to be defined as 0 to 360. In other words Greenwich would be at 180 degrees longitude.


texts: Controlling "Fill Properties" in ArcGIS, so that that symbols are visible in legend?

texts:
I use marker fill symbols in my legend in ArcView 10, in order to identify different habitats, and frequently have to increase the separation of the symbols in order to make the map more readable.  However, increasing the separation causes the legend to sometimes not include the marker symbol (ie it falls off the edge of the legend patch - see attached where "SI" is only partially visible).

I was wondering if people had come across this problem and new of a solution, I have used the following but find them to be a temporary rather than permanent solution.

(A) Increase the Patch Size - not always possible when there are lots of labels and not much map space.

(B) Use "Symbol Property Editor" -> "Fill Properties" and change the "Offset" x and y values so that the marker symbol appears in the legend.

(C) Use "Convert to Graphics", this does not work because when you move the label the marker disappears and also you lose the interactive functionality of the legend.

(D) Use another type of symbology, but this is problematic because I am following a industry standard palette and no colour scheme would work with the large number of categories that I have.


texts: arcpy.da.SearchCursor, cannot read preceding zeros from csv

texts:
I'm using arcpy.da.SearchCursor to read a csv that has a column with preceding zeros, e.g.:

myCSVColumn
0124578
12598
04589
012568
...

When using the arcpy.da.SearchCursor it returns the value w/o the preceding zeros:

124578
12598
4589
12568

In addition, when using Table to Table method it returns the same result.  Does anyone know a workaround for this besides opening the csv with python and parsing through it to get the values?


texts: ogr2ogr + MapInfo TAB + SQLite (SpatiaLite) = unicode problems

texts:
I'm adding MapInfo TAB files to SQLite (SpatiaLite) database using those arguments:

ogr2ogr -f "SQLite" dataset.sqlite somelayer.tab -dsco spatialite=yes

Source file somelayer.tab have WindowsCyrillic encoding defined:

!table
!version 450
!charset WindowsCyrillic

Definition Table
  Type NATIVE Charset "WindowsCyrillic"
  Fields 1
    Name Char (100) ;

As a result, dataset.sqlite have somelayer table with VARCHAR Name column.

But I need store Name values in unicode.

How can I tell ogr2ogr use unicode NVARCHAR type for Name column and encode text from WindowsCyrillic (Windows-1251) to UTF-8 encoding for storing.


texts: Fiona error adding property to shapefile

texts:
I'm getting the following error when trying to add a property to a .shp file using fiona in python:
    Traceback (most recent call last):

  File "<ipython-input-38-5f04c3b56601>", line 17, in <module>
    j=j+1

  File "/Users/alex/anaconda/lib/python2.7/site-packages/fiona/collection.py", line 402, in __exit__
    self.close()

  File "/Users/alex/anaconda/lib/python2.7/site-packages/fiona/collection.py", line 386, in close
    self.flush()

  File "/Users/alex/anaconda/lib/python2.7/site-packages/fiona/collection.py", line 376, in flush
    self.session.sync(self)

  File "fiona/ogrext.pyx", line 939, in fiona.ogrext.WritingSession.sync (fiona/ogrext.c:15649)

RuntimeError: Failed to sync to disk

Here is my code - I'm adding a dummy property 'jazz' for now but I have many properties with both float and string types that I want to add once I figure out how to get it to work...

import random
#
import fiona

with fiona.open("data/subplace_maps/SP_SA_2011.shp") as input:
    schema2 = input.schema.copy()
    schema2['properties']['jazz'.decode('utf-8')] = 'str:50'
    temp = input.crs
    with fiona.open('data/subplace_maps/output.shp', 'w', crs=input.crs, driver=input.driver, schema=schema2) as sink:
        # iterate over shapefile shapes
        j = 0
        for elem in input:
            # update progress
            if(j % 500 == 0):
                print j,'/',len(input)
            #
            ja = ('jazz' + str(random.random())).decode('utf-8')
            elem['properties'].update(jazz=ja)
            #
            sink.write(elem)
            j=j+1

print 'done'


texts: Point and Polygon Data refusing to display: corrupt data?

texts:
I have two sets of shapefiles, one point data and one polygon data. Both were sent as part of a bigger set of data.

Neither the point or polygon data will display when loaded into QGIS. The raw data has not been edited, the CRS is correct for the project and layer (27700). I have tried multiple times to load these layers with no success. I have have remaped the Spatial CRS from within the layer properties.

I am at my wits end trying to work out what is wrong!


texts: can geoserver select one of many rasters from a postgis table

texts: I have a postgis table with many different rasters in it. Is it possible for geoserver to grab only one of the rasters from the table based on something like a where clause?


texts: Missing shp2pgsql command on postgis docker installation, when I use postgis:postgis:12-3.1 image?

texts:
I'm installing a docker container from this image, taken on Docker Hub : postgis:postgis:12-3.1, and according to posts on Stackexchange I did a service reload after the container initialization.
But the shp2pgsql of Postgis isn't here. A locate cannot find it, inside the container.
Must I do additional steps?


texts: Arcgis python： parameter not changed leads to a dead loop

texts:
Requirement Analysis：
 I have got a Tin/Terrain file, and I set an initial water level, the target is to calculate the volume of this place under the initial water level using the Terrain file. Comparing the volume with the actual value, I got the ideal water level.

Tools: Arcgis 10.1 (polygon volume_3D module) python 2.7

code:

import arcpy
import os

arcpy.env.workspace = "G:/0arcpy/hjh_split"
arcpy.env.overwriteOutput = True
arcpy.CheckOutExtension("3D")

fclist = arcpy.ListFeatureClasses()

inSurface = "G:/0arcpy/hjh_tin"
zField = "sswl"
refPlane = "BELOW"
sAreaFld = "SArea"
volFld = "Volume"

for fc in fclist:
    arcpy.PolygonVolume_3d(inSurface, fc, zField, refPlane, volFld, sAreaFld)
    rows = arcpy.UpdateCursor(fc)
    for row in rows:
        i = row.getValue("Volume")
        j = row.jyvol
        while ((abs(i-j))/j>0.1) :
            row.setValue("sswl",row.getValue("sswl")+0.01)
            rows.updateRow(row)
            arcpy.PolygonVolume_3d(inSurface, fc, zField, refPlane, volFld, sAreaFld)
        print "i="+str(i)+"-"+"sswl="+str(row.sswl)
        print "Finished"

Result:

Question:
 The water level(sswl) is changed with plus 0.01, but The i(volume) value is not changed leading to dead loop. But why the value i = row.getValue("Volume") is not changed in the loop. What's the problem and how to change the code?


texts: Creating shapefile in standalone script using PyQGIS?

texts:
I'm trying to write a standalone script that writes a shapefile (using QGIS 3.0.2). I've tried something like :

from qgis.core import *
from qgis.PyQt.QtCore import QVariant
from qgis.utils import QGis

QgsApplication.setPrefixPath("/usr/bin/qgis", True)
qgs = QgsApplication([], False)

qgs.initQgis()

fields = QgsFields()
fields.append(QgsField("ID_HYD", QVariant.Int))
fields.append(QgsField("SBD_KM2", QVariant.Double))

writer = QgsVectorFileWriter("/home/sylvain/test.shp",
                             "CP1250",
                             fields,
                             QGis.WKBPoint,
                             None,
                             "ESRI Shapefile")

qgs.exitQgis()

The issue is as follow :

/usr/bin/python3.5 /home/sylvain/test.py
Must construct a QGuiApplication first.

Process finished with exit code 1

Is there a way to write a shapefile on drive without constructing a GUIApplication ?

If not, how should I do do create such an item ?


texts: How to know what is the correct ESPG code for clipping Landsat imagery?

texts:
This is somewhat related to a previous question of mine.

My goal here is to clip Landsat imagery to a rectangle that I've defined by four points of latitude and longitude in WGS84. To do this, I've used this python module to convert the points of the rectangle into utm coordinates suitable for use as arguments to gdalwarp. My question is how do I determine what the correct EPSG code is for gdalwarp? I assumed it would be a simply mapping of UTM Zone to EPSG (for example, coordinates in UTM Zone 11N (southern california) would map to http://www.spatialreference.org/ref/epsg/wgs-84-utm-zone-11n/). However, when I try to use this code in my gdalwarp command, it hangs trying to create an image much too large. When I change the EPSG code to http://spatialreference.org/ref/epsg/nad83-utm-zone-11n/ it works and gives me the clipped imagery I expect. What I don't understand is that the info page on the USGS web site implies that 32611 is the correct EPSG code for this specific landsat tile, but it seems that 26911 is the actual correct EPSG code. And further, why don't either of these EPSG codes appear in the output from gdalinfo for the file?

The command I'm using to clip the images is:

gdalwarp -of gtiff -t_srs EPSG:<EPSG Code goes here> -te <bounds in utm>  input.tif output.tif

The output from gdalinfo for my input image is:

Size is 15381, 15681
Coordinate System is:
PROJCS["WGS 84 / Pseudo-Mercator",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0],
        UNIT["degree",0.0174532925199433],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Mercator_1SP"],
    PARAMETER["central_meridian",0],
    PARAMETER["scale_factor",1],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    EXTENSION["PROJ4","+proj=merc +a=6378137 +b=6378137 +lat_ts=0.0 +lon_0=0.0 +x_0=0.0 +y_0=0 +k=1.0 +units=m +nadgrids=@null +wktext  +no_defs"],
    AUTHORITY["EPSG","3857"]]
Origin = (-13185978.181983433663845,4059663.631883424706757)
Pixel Size = (18.129383469292915,-18.023985075689048)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (-13185978.182, 4059663.632) (118d27' 5.97"W, 34d13'46.30"N)
Lower Left  (-13185978.182, 3777029.522) (118d27' 5.97"W, 32d 6'16.03"N)
Upper Right (-12907130.135, 4059663.632) (115d56'48.20"W, 34d13'46.30"N)
Lower Right (-12907130.135, 3777029.522) (115d56'48.20"W, 32d 6'16.03"N)
Center      (-13046554.158, 3918346.577) (117d11'57.08"W, 33d10'24.35"N)
Band 1 Block=15381x1 Type=Byte, ColorInterp=Red
  NoData Value=0
Band 2 Block=15381x1 Type=Byte, ColorInterp=Green
  NoData Value=0
Band 3 Block=15381x1 Type=Byte, ColorInterp=Blue
  NoData Value=0

Sorry if this is a dumb or obvious question, I'm still new to GIS.

[EDIT]

Finally figured it out!

I was using landsat-utils to process the landsat imagery down to a single rgb tiff. In the process, the image was getting mapped to web mercator from the UTM of the original input imagery. I'm still not sure why I was getting reasonable results with UTM, but by using this gdalwarp command, everything comes out more or less as expected.

gdalwarp -of gtiff -t_srs EPSG:3857 -te <bounds in web mercator>  input.tif output.tif


texts: MultiPolygon Intersect only working for first Polygon in collection, how can I fix this?

texts:
I am very new to PostGIS.

I have a table called regions, that has name, an id and a geojson that is a MultiPolygon.

Now I need to check whether certain points or other MultiPolygon regions intersect with it.

This is my sql statement:

SELECT regions.id, name
        FROM regions
        WHERE ST_INTERSECTS(
          regions.geom,
            ST_SetSRID(ST_Makepoint(lat,long), 3857));

This works fine for the first Polygon of the MultiPolygon of a region. But if I check if a point/MP intersects with any of the other MultiPolygons within a region it returns nothing.

Could anyone tell me how to fix this?


texts:
I'm using snapping in a web map for a move/rotate custom tool. According to the API documentation the default option is to set all feature and graphics layers in the map to be target snapping layers.

The snapping seems to work fine since the user is getting the X symbol over the vertices of the features when the pointer is close with the CTRL key pressed.

However when I capture the point coordinates on click, the value is not identical to the feature's vertex coordinate as I expected.  There's a small difference and some variance in multiple clicks for the same vertex.  For the tool I really need to have the same exact coordinate of the feature so the operation will no generate gaps in the layer.

I'm using ArcGIS API for JavaScript 3.7 and the drawEnd event (using esri.toolbars.Draw.POINT) to get the coordinate.

Is this the expected behaviour for the snapping (an approximation instead of the exact coordinate of the feature) ? Or is there a way to get what I want.


texts:
I'm trying to find the best route through multiple points, considering traffic, speed limit, one way, and so on. ArcGIS online has Plan Routes functions which we only need to provide point feature class to get best routes that consider traffic and other variables. However, ArcGIS Pro does not provide that function. We have to use Network Analysis where we have to create network dataset first to have best routes that consider traffic and speed limits.

Is my understanding right ?


texts: Showing multiple burials in same grave for Cemetery GIS?

texts:
I am an intern creating a Geodatabase for a municipal cemetery.

The one lot can hold up to 3 people (2 people, 1 cremation)  two people at the most.

How would I show that information on the map using a polygon while including the position of each body?

The lots are a fixed size. I am using ArcMap 10.2.2


texts: Can I transfer catagorised symbol groupings to to an attribute?

texts:
I have a layer file in ArcGIS 10.2 that symbolizes some points by category. 100s of different values in the value field have been grouped down to about 20 symbols. Unfortunately the data I want is these symbols. I can see this but I want to be able to read it as an attribute or select all points that share a symbol?

Is there anyway I can export the symbol definition to a table or select according to symbol?


texts: What driver should I use when trying to open a Geospatial PDF with OSGEO and Python?

texts:
My code is as follows:

from osgeo import ogr
from osgeo import gdal

dso = ogr.Open('myfile.pdf')
print dso

dsg = gdal.Open('myfile.pdf')
print dsg

driverGeoPDFogr = ogr.GetDriverByName('PDF')
dso2 = driverGeoPDFogr.Open('myfile.pdf')
print dso2

This results in:

None
ERROR 4: `myfile.pdf' not recognised as a supported file format.

None
None

I know the file is a valid file, and when running gdalinfo from command line, the result for driver is:

Driver: PDF/Geospatial PDF

My question is, is PDF the incorrect driver name here for Geospatial PDF, or is this file format incompatible here?  I'm running Python 2.6 and GDAL 1.11.2.  Note, this is a Geospatial PDF, NOT a GeoPDF.


texts: Generate Cesium 3-D tiles from shapefile with FME

texts:
I'm trying to generate 3-D tiles from a shapefile. Just simple shp reader directly to Cesium 3D tiles writer.  It does generate files,  but after uploading them to Cesium Ion it does not work.

It can't even generate preview of it.

Am I missing some step in FME?

Additional info:

While adding SHP reader I used "Read from source" option, that's why it is blank now:

while adding 3-D tiles writer I used "Same as source" option.

Settings of shapefile:

I also tried to use EPSG:4326 in both components but same result.


texts: Determine 'native' size of wms layer based on GetCapabilities response

texts:
I'm trying to understand the relation between a <LatLongBoundingBox> and <BoundingBox> in a GetCapabilities WMS response, given that the SRS of the bounding box is EPSG:4326.

Here's an example of the WMS that I'm trying to query:

<Layer>
<!-- snip -->
<LatLonBoundingBox maxx="37.968" maxy="70.144" minx="-4.646" miny="48.833"/>
<BoundingBox SRS="EPSG:4326" maxx="0.662672" maxy="1.224240" minx="-0.081091" miny="0.852295"/>
<!-- other BoundingBoxes omitted -->
</Layer>

My questions related to this response:

Why are the BoundingBox values different from the LatLonBoundingBox values an how is this useful? It looks like a scaling hint of some sorts.
Is it possible to determine the 'native' resolution of the layer (in pixels) using these values?


texts: QGIS plugin external libraries in plugin’s directory

texts:
I’m developing my first QGIS plugin and I have to use two external libraries not included in the default OSGeoW64 installation: Couchdb and Pandas.

I would like to include those libraries in my plugin’s directory so that users don’t have to deal with the installation of them.

So I created a folder called lib, inserted in it an empty _init_.py and downloaded in it the libraries using pip:

python3 -m pip install -t lib couchdb

The directory looks like this:

mhyconos
|--- lib
     |--- bin
     |--- couchdb
     |--- CouchDB-1.2.dist-info
     |--- _init_.py
|--- _init_.py
|--- mhyconos.py
...

In mhyconos.py I imported couchdb using

from lib import couchdb

but when I install the plugin, this is what happens:

What am I doing wrong?


texts: qgis-bin-exe-stopped-working No Shapefiles allowed

texts:
I'm using Windows 7 and I'm having the problem of QGIS 2.18.16 crashing with any shapefile. Rasters seem to be OK. I have tried uninstalling and reinstalling 2.18.15, 2.18.14 and one older version (2.17 or earlier, I can't remember). I turned off all plugins before trying anything. I found earlier advice about changing a registry key, so I used the Windows command line to start regedit and under HKEY_CURRENT_USER> SOFTWARE> QGISI deleted all the QGIS and then did the 2.18.16 install. I also downloaded the OSGeo4W version and tried running QGIS from that.

I haven't tried any other vector formats because I really need to work with shapefiles.  I'm not having any trouble with my desktop, just my laptop. I was using the style editor when the trouble first appeared.

I looked all over at the other posts about this problem, like qgis-bin.exe stopped working and I'm not sure how to provide more information.

Edit: I tried loading a geodatabase with vectors, it worked and I changed a few symbols and saved. Upon exiting I received this same message and now I get this fail when I try opening geodatabases. I also get this fail when I'm loading the new project (with just this simple geodatabase loaded).


texts: Trapping error messages with arcpy.getmessage issue

texts:
I am trying to capture an error with arcpy.getmessage(2) within a try/except statement but when it fails it tells me the issue is with the getmessage, it is not returning the error that occurred in the try part. This is what I have setup to capture the error message...

try:
    arcpy.CopyFeatures_management(Flow_Arrow, Flow_Arrow__2_, "", "0", "0", "0")
    print "Flow Arrow Updated"
    print "Flow Direct Feature Set Updated"
except:
    print "Error"
    print arcpy.getmessage(2)

I know it has got to be something simple I am just overlooking, I tried arcpy.getmessage(), and arcpy.getmessage(0), I went with the (2) because I was googling the issue it was one of the proposed solutions.


texts: Read BigGeoTIFF multiband GDAL using Python

texts:
I'm trying to read\load BigGeoTIFF multi-band raster (>4GB), using GDAL module in python but I can't seem to make it. At first the error was that BigGeoTIFF format is not supported by GDAL and libtiff versions, but according to this document I have the required versions. I also saw somewhere (can't find it now for some reason) that you need to specifically enable the BigGeoTIFF format in GDAL but I don't know how to do so. Anyways, now when I try to read the raster with:

source = r"C:\Users\Administrator\Desktop\california_merge.tif"
dataset = gdal.Open(source)
print type(dataset)
>>>NoneType

When I use this code on <4GB multi-band raster I get the desired result

source = r"C:\Users\Administrator\Desktop\california_merge_clip.tif"
dataset = gdal.Open(source)
print type(dataset)
>>>osgeo.gdal.Dataset

What can I do to be able to read (and later process) BigGeoTIFF file in Python (version 2.7 on win 10)


texts:
I want to create Subtypes for a feature class, but the option is not active.
Under Feature Class Property I created some Fields, but the Subtypes tab is not active.
I can't choose Subtype Field and can't type under Code or Description.
I did try to create subtypes for another map, even open a new document and created a new File Geodatabase, new feature class, but Subtype Tab is not working.
I am using ArcMap 10.7.


texts:
Is it possible to load an .las file directly into QGIS?

I would like to be able to load LiDAR data and to be able to do analysis based various feature attributes; e.g., elevation, classifications, etc. While LASTools (which was suggested in What is the procedure to load LAS files in QGIS 2.0.1?) is a nice set of tools, it still didn't produce an .las point file that could be read by QGIS


texts:
I'm adding a label to a Print Layout using python. I do so with the following code:

label = QgsLayoutItemLabel(layout)
label.setText("This is my map label")
label.adjustSizeToText()
layout.addLayoutItem(label)
layout.attemptMove(QgsLayoutPoint(0.25, 0.25, QgsUnitTypes.LayoutInches))   #repositions map label

The label is small font by default and I want to resize it, but cannot figure it out. After looking at the documentation for QgsLayoutItemLabel() there is no method to set the text label size.

The problem here is label.setText(), creates a string with the text.  I could modify the size of the label but only if the string is some kind of text object?

In this example in the first answer, the person created a QgsPalLayerSettings() object, which he then formatted.

I found a helpful example in the source code for the QuickPrint plugin. The author gets user input via a dialog box using the following line:

titel = self.dlg.titelFld.text()

Here he converts this input to text using the text() method, but I can't figure out how to apply text() method to a string.

Any suggestions?


texts:
I would like to work with the USGS dataset PADUS in R. However the files are huge and I just want to plot out the main state and national forests, parks and protected areas. Does anyone have any suggestion on how to load and plot these very large files?
USGS DATA DOWNLOAD

Normally I would use readOGR (rgdal) to load in the files but this just seems to crash or fail to load.


texts: Installing QGIS on computer which has ArcGIS Desktop?

texts:
I have ArcGIS Desktop on my computer, but I want to install QGIS 3.8 as well. The installation finishes successfully as it says, but I can't find many files after and the program is not working.
I think, the problem is ArcGIS Desktop 10 which I have on the same comuter. I tried to install QGIS 3.2 too, with the same bad result.

My solution was that I had a clean QGIS install on another computer and copied those files to this installation folder, and the program is now working.

Is there any interference between these (ArcGIS Desktop & QGIS) programs?

Even with Python or something else?


texts: How to substring a number out of a string field using Python field calculator in ArcGIS?

texts:
Using ArcGIS Desktop 10.7
I want to substring numbers as text out of a string field into a new string field with field calculator so that I can use them for layer symbology. The field contains highway names such as 'I- 280' and 'State Rte 25' and I want only the '280' and '25' for a new field. Is there a Python conditional statement I can use that will substring the number regardless of digit length, placement (end or between text, eg, "State Rte 25 Bus")?


texts: Disperse markers - trouble with input file?

texts: I had an excel sheet listing addresses that I imported into ArcGIS and geocoded which created a File Geodatabase Feature Class of all the points. I am trying to display the map where it disperses points that are overlapping, but it keeps telling me the input is wrong for the tool Disperse Markers (Geocoding_result.lyr does not contain representations)...to what and how do I convert the geocoded results to have the correct input to use for the tool Disperse Markers? Using ArcGIS 10.5


texts:
Dear all am trying to perform  a union between two shapefiles ( for the moment ) - line 88-

 oLayer.Union(oLayer2,oLayer3)

and I'm getting stack with this error

Traceback (most recent call last):
  File "/home/Uka/PycharmProjects/histcult/histoire.py", line 88, in <module>
    oLayer.Union(oLayer2,oLayer3)
  File "/usr/lib/python3/dist-packages/osgeo/ogr.py", line 2777, in Union
    return _ogr.Layer_Union(self, *args, **kwargs)
RuntimeError: Attempt to write non-polygon (LINESTRING) geometry to POLYGON type shapefile.

and this is my code :

import os
from osgeo import ogr
import osgeo.osr as osr

ogr.UseExceptions()
osr.UseExceptions()

# setting up the shapefile  driver
shpDriver = ogr.GetDriverByName("ESRI Shapefile")

# creating the data source
OutputDS = shpDriver.CreateDataSource("test_sample/output_histoire.shp")
OutputDS2 = shpDriver.CreateDataSource("test_sample/output_histoire2.shp")
OutputDS3 = shpDriver.CreateDataSource("test_sample/output.shp")

# create the spatial reference, WGS84
srs = osr.SpatialReference()
srs.ImportFromEPSG(2154)

# creating output layer the layer
oLayer = OutputDS.CreateLayer("culture", srs, ogr.wkbPolygon)
oLayer2 = OutputDS.CreateLayer("culture2", srs, ogr.wkbPolygon)
oLayer3 = OutputDS.CreateLayer("output", srs, ogr.wkbPolygon)

#Output Field names
cultureField = ogr.FieldDefn("cultPr", ogr.OFTString)
cultureField.SetWidth(30)
oLayer.CreateField(cultureField)
oLayer.CreateField(ogr.FieldDefn("ID_C", ogr.OFTInteger))

#input layers
cultureUrl = "test_sample/culture.shp"
parcelleUrl = "test_sample/parcelles2016test.shp"

cultureDtSource = shpDriver.Open(cultureUrl, 0) # 0 means read-only. 1 means writeable.
parceDtSource = shpDriver.Open(parcelleUrl, 0) #
cultureLayer = cultureDtSource.GetLayer(0)
parcelleLayer = parceDtSource.GetLayer(0)
parcelleDef = parcelleLayer.GetLayerDefn()

# Check to see if shapefile is found.
if cultureDtSource is None  :
    print ('Could not open %s' % (cultureUrl))
elif parceDtSource is None :
    print('Could not open %s' % (parcelleUrl))
else:
    print ('Opened %s' % (cultureUrl))
    print('Opened %s' % (parcelleUrl))

    cultureFC = cultureDtSource.GetLayer().GetFeatureCount()
    parcelleFC = parceDtSource.GetLayer().GetFeatureCount()
    print ("Number of features in %s: %d" % (os.path.basename(cultureUrl), cultureFC))
    print("Number of features in %s: %d" % (os.path.basename(parcelleUrl), parcelleFC))
    for i in range(parcelleDef.GetFieldCount()): # iterate on parcelle features to get the names
        print (parcelleDef.GetFieldDefn(i).GetName())
    nbrIntersection = 0

    for featureP in parcelleLayer: # iterate on parcelle features to get the the geometry

        geomP = featureP.GetGeometryRef()
        for featureC in cultureLayer: # iterate on culture  features to get check the intersection
            geomC = featureC.GetGeometryRef()

            intersection = geomC.Intersection(geomP)
            if not intersection.IsEmpty():
                # create the feature
                outputFeature = ogr.Feature(oLayer.GetLayerDefn())
                # Set the attributes
                outputFeature.SetField("cultPr", featureC.GetField("cult"))
                # Set the feature geometry
                print(intersection.GetGeometryName())
                outputFeature.SetGeometry(intersection)
                # Create the feature in the layer (shapefile)
                oLayer.CreateFeature(outputFeature)
                # Dereference the feature
                feature = None
                nbrIntersection = nbrIntersection + 1
        cultureLayer.ResetReading()
    parcelleLayer.ResetReading()
    parcelleLayer.Erase(oLayer,oLayer2);
    # add union between oLayer and oLayer2
    oLayer.Union(oLayer2,oLayer3)
    print('you got ', nbrIntersection, ' intersections')
    # Save and close the data source
    OutputDS = None
    OutputDS2 = None
    OutputDS3 = None


texts: How can GIS be used in disaster management?

texts:
I know that mapping plays a significant role in dealing with natural disasters, but I'm not familiar with how that works as my background is primarily web development.

I'm wondering if there is any way for GIS professionals or amateurs to contribute to the relief efforts now or as part of recovery efforts.


texts: How to draw a Circle of fixed radius in QGIS using QgsMapCanvasItem?

texts:
I have derived a class from QgsMapCanvasItem for drawing circle of some radius at a given latitude and longitude.
In my paint function I am using QPainter to draw the circle as

void paint (QPainter* painter){
 painter->drawEllipse(center, radius, radius);
}

void updatePosition(){
   center = toCanvasCoordinates(center_lat_lon);
}

The circle is drawing properly and its position is getting updated whenever the events like panning, zoom in, zoom out are happening.
My problem is, whenever I am zooming in or zomming out, radius is not changing accordingly.  Also I want to update my radius according to the projection system  ( say WGS84 ). Any help is greatly appreciated.


texts:
I have two tables with linestrings, the first one represent wires and the second the infrastructure containing these wires.

While drawing the cables, it can happen that the wire's vertices do not match perfectly the infrastructure ones, it is generally a matter of a few cm.

Is there a way to update the whole wire table geometries by making that each vertex is matching the closest one from the infrastructure geometry?

EDIT :

Here is a picture to show the situation. In red is my infrastructure, in blue the wire. As we can see on the upper case the middle point of the wire is not perfectly matching the one of the infrastructure, and it should.

The case here is simple, in reality I have full networks of each. I would like to automatically update the wire geometries so that they match the infrastructures. Here I would like to find a query which would do like in the picture.

I tried ST_Snap, but I cannot manage to make it work. I tried something like that, which didn't give any conclusive result :

UPDATE wire
SET geom = ST_Snap(wire.geom, infra.geom, ST_Distance(wire.geom, infra.geom)*1.5)
FROM infra

EDIT 2:

I must precise that unlike in my example, a wire can be passing over more than one infrastructure, so if there was a way by replacing the wire's geometry by the infrastructure one, it wouldn't work properly in my case.


texts:
I am working with this query:

SELECT ST_MakeLine(route.geom ORDER BY seq)  as test
    FROM (SELECT * FROM pgr_fromAtoB('way',48.24287652969359,33.85220534163203,48.24991464614867,33.85092228417727  )) As route

the above query returns a null value. I tried:

SELECT * FROM pgr_fromAtoB('way',48.24287652969359,33.85220534163203,48.24991464614867,33.85092228417727  )

this query works fine and returns a set of columns with a geom (line) and seq columns witch is sequence of parts of lines.

I think I am doing something wrong with ST_MakeLine but can not find where it is.
Can you please help me find the mistake? Does St_makeline work fine with a series of line geometries?


texts: Having trouble getting geojson to show up with leaflet from geoserver

texts:
I have enabled jsonp in geoserver/WEB-INF/web.xml. It appears in my layer preview in geoserver. I cannot seem to get my geojson layer to show up in my map. The tiles show up but with no layer. Not sure where to go from here. I have been searching through other peoples coding but no luck. Wondering if anyone could see my problem here with my code.

  <title>Geojson Leaflet</title>
  <link rel="stylesheet" href="http://cdn.leafletjs.com/leaflet-0.7.3/leaflet.css"/>


texts:
I have a visualization in CartoDB, and I've added a Title Item with a url with the format [Name] (url). This url opens in the window of cartodb.

Is it possible to open the url in a new window?

Is it possible to add charts in cartoDB?

I have heard about Highcharts or Googlecharts, but I want to add the chart in the window of cartodb, I don´t have a own web page of html.


texts: Building footprint polygon from LoD1 CityGML 2.0 using R and sf?

texts:
I have LoD1 building data from a federal agency in Germany in the CityGML format. I would like to import and strip the data into R to get only the footprint of the buildings and an attribut about the height. In case of building parts with different heights I would ideally get a the building footprint with multiple height fields (or a list containing the different height information). sf let's me import the data but I only partially gives the expected result: While it works for some buildings, for other building only a part is imported. I included a sample from the data with one building where the import works just fine and the other where only parts of building are imported. I noticed that the structure of the cityObjectMember from the first bldg:Building is different and the gml:id for this building is empty.  Although I guess this is the reason why it's not imported correctly but I was wondering if this is because the data is corrupted or if there is a logic behind this.

I have no experience with CityGML but I also couldn't find tutorials on how to work with CityGML in R or any other editor that would allow for anything beyon visualization.

How do I import and strip CityGML data in R?

Here's what I've tried so far:

Test data with two buildings

library(dplyr)
library(sf)

#Import buildings
buildings = read_sf("test_building.xml")

# First building (with empty gml_id) not correctly imported
# Only spatial information for a part of the building is imported
plot(st_geometry(building[1,]))

#Second building (with gml_id==DENILD4300007ABe) works fine
plot(st_geometry(building[2,]))


texts: MMQGIS Geocoder Request Denied Error - Using Google API Key

texts:
I need to geocode 12,000 addresses, so I enabled billing and set up credentials to use the Google Geocoding API.  When I try to run the MMQGIS geocoder using my new API Key, I immediately get the message: Request Denied - 0 addresses geocoded.  I have enabled the Google Maps Geocoding API as well as the Google Maps Javascript API.

What else do I need to do to get this working?


texts:
I'm aiming to represent a huge geojson file (around 250Mb) using leaflet, therefore, I need to render my data before displaying it.
I want to only load data that is in the bounds of my leaflet frame. Also, I want, as long as the map load the quality of the polygon bounderies increases. I think the same process is used in loading tiles.


texts:
I have QGIS 2.4 installed on a machine running Windows 7 64-bit "by accident". I was using 2.2 and when I added some package(s) it also installed the 2.4 version.

Can anyone provide some advice with the following?

Remove version 2.4 and ensure there are no hooks, etc remaining on my computer (e.g. best method such as using Revo uninstaller or other)
Re-install version 2.2 (I'm assuming it would be download it from here:
http://qgis.org/downloads/ , specifically the file
QGIS-OSGeo4W-2.2.0-1-Setup-x86_64.exe, dated February 26, 2014)
Ensure that adding packages in the future doesn't automatically upgrade
my version 2.2 to 2.X.


texts: Computing the density for each layer with lidR

texts:
I am trying to perform a 3D profile index as a function. For this, I need to slice the layers and calculate the point density for each sliced layer.  I write a for loop;
The problem is
density=grid_density(subset,res=0.1)

and
pt<-as.raster(pt)

parts gives always an error. When I run  as.raster I always get this error:
Error in UseMethod("as.raster") :
  no applicable method for 'as.raster' (applied to an object of class "c ('LAS', 'Spatial')")

I tried raster() instead of as.raster() it works but I don't know what is the main difference.
I filter the points with filter_poi() . In this filter I specified z value, the points must be calculated lower than the next sliced layer and higher than the layer itself.  Then I want to calculate the grid_density() it gives this error:
Error in validityMethod (object): invalid extent: xmin> = xmax
Additionally: Warning messages:
1: In min (x @ data $ X): No missing arguments for min; Returning inf
2: In max (x @ data $ X):
   No complete arguments for max; Returning inf
3: In min (x @ data $ Y): No missing arguments for min; Returning inf
4: In max (x @ data $ Y):
   No missing arguments for max; Returning inf

library(lidR)
library(raster)
library(sp)
library(rgdal)
library(gridExtra)

#set the parameters
sz=0.05 #slice size in meters
s=seq(0,2,sz) #number of layer
k=seq(-1, 2.00,0.05) #correction factor (k)

las <- readLAS(x)
pt=grid_density(las,res=0.1) #pt is the total LiDAR points for all layers
pt<-as.raster(pt) #convert to raster
crs(pt) <- "+proj=utm +zone=31 +datum=WGS84 +units=m +no_defs"

# #make base layer, used to maintain the extent of the Pi raster layers
base_layer<-pt
values(base_layer)<-0

density_rasters=c()
for (i in s){
  subset=filter_poi(las, Z>=i & Z<(i+sz))
  if (!is.null(subset)){
    density=grid_density(subset,res=0.1)
    ras=as.raster(density)
    setExtent(ras, ext=extent(base_layer),keepres = TRUE, snap = TRUE) # set extent to baselayer, otherwise the extent will be to small to be able to stack the rasters
    ras=merge(ras,base_layer)
    density_rasters=c(density_rasters,ras)

  }
  else { #used to insert empty raster if no points are within a given pi layer
    density_rasters=c(density_rasters,base_layer)

  }
  density=0
  ras=0
  subset=0
}


texts: Testing for spatial correlation using two point layers in QGIS?

texts:
I have two gridcell-based point layers that I am interested in investigating their spatial correlation in QGIS, or alternatively in R as well. The variables I am interested in are:

Map 1: Point layer of percent change in price
Map 2: Point layer of demand estimate

A few word on the context of the analysis. The maps are generated via a simple supply-demand pricing framework, which is spatially explicit at the gridcell level. For each gridcell, we have a supply curve that characterizes the availability of a resource k based on data on availability and procurement cost. Using scenarios from a GIS-based model, I can generate different demand schedules, of which Map 2 is an example. First, we use a business-as-usual scenario to calibrate a BAU price vector. Then, changing certain constraint parameters in the GIS-based model, we run alternative scenarios which generate a different demand profile, and for which we generate their specific price vectors using the supply curve. Thus, the price impacts represent the percent change deviation from the BAU and the alternative scenario of which Map 1 is an example.

Hence, what I am interested in is to somehow evaluate the degree of correlation between the demand location and the price impact location. My questions are:
What is the best way to go about this?
Is it possible to do correlation analysis in QGIS?
Or is best to run spatial regressions to test the fit between the price impact and demand layers?

I have looked at certain alternatives in QGIS. Specifically, I have managed to perform a Nearest Neighbour Analysis (NNA) to investigate the linkages between the location of the price impacts in Map 1 and the location of the demand in Map 2. However, I wonder if that is the 'best' route to go about.

Any advice and/or pointers would be very welcome.


texts: Generating tiles for osm/esri shapefiles

texts:
I've generated tiles for OpenStreetMap .OSM files using Maperitive and quickly consumed them with leaflet (with no tile cache yet).
I'm now looking for the right and simple way to do that also with the shapefiles given by OSM.

(Tried TileMill but the generated mbtiles file needs further processing)


texts: tiny JS discretization library for choropleth representation

texts:
I would like to know if a tiny stat/discretization js library ready to use exists.

I do know about the existence of jStat (which I may explore a little bit more thoroughly), and I know that mapfish has a geostats module/widget to create choropleth map - I found the code here. But it relies on openlayers and I would like to have a standalone lib and eventually be free to communicate with openlayers.

I would like to have few basics methods:

min
max
mean
median

and be able to implement the following discretization with the ability to specify the number of classes and the precision (decimal):

Quantile
Equal Interval
Natural Breaks (Jenks)

If it does not exist yet, I will probably code something myself and release it on github under free license. Any resource and advice is welcome.


texts: Setting attribute values using PyQGIS in an edit session

texts:
My PyQGIS script works fine if I am not in an edit session, but 99% of the time I will be, so I need it to work whilst in editing mode.
from PyQt5.QtGui import *

layer = QgsProject.instance().mapLayersByName("test")[0]
iface.setActiveLayer(layer)

value='222'
with edit(layer):
    for feature in layer.selectedFeatures():
        feature['flag1'] = value
        layer.updateFeature(feature)

msg = QMessageBox()
msg.setText("All done.")
msg.show()

So basically, it just updates the 'flag1' field with the assigned value for those features that have already been selected from the 'test' layer. How do I make this work inside an edit session? Currently I receive this error:

File
"C:/PROGRA~1/QGIS3~1.16/apps/qgis/./python\qgis\core\additions\edit.py",
line 38, in enter
assert self.layer.startEditing() AssertionError


texts:
Below is so far I've done to convert .csv file into shapefile

Mydata <- read.csv("londonburglary201402.csv", header = T, sep = ",")
class(Mydata)
coordinates(Mydata) <- ~Longitude + Latitude
class(Mydata)
[1] "SpatialPointsDataFrame"
attr(,"package")
[1] "sp"

but I don't know how to write the shapefile with projection.
I am trying to use rgdal packages to do.


texts:
I need to convert a bunch of coordinates that were stored in Feet to Decimal Degrees. I found a Function in VB.Net, though it's results are not coming back correctly. When I input 3106972, 1417936 into this function I get: -80.7964701731349, 32.8465039721407 where they should be: -105.115621, 40.480104. Where am I going wrong here?

Here is the function:

Private Function ConvertFeetToDecimalDegrees(lat As Integer, lng As Integer) As String
' This function converts coordinates from Florida State Plane North Zone NAD83 feet
' to Geographic coordinates in decimal degrees
' Set up the coordinate system parameters.
Dim A = 20925604.47                         ' major radius of GRS 1980 ellipsoid, feet
Dim Ec = 0.0818191910435                    ' eccentricity of GRD 1980 ellipsoid
Dim Ec2 = Ec * Ec                           ' eccentricity squared
Dim AngRad = 0.0174532925199433             ' number of radians in a degree
Dim Pi4 = 3.141592653582 / 4                ' Pi / 4
Dim P1 = 29.583333333333329 * AngRad        ' latitude of first standard parallel
Dim P2 = 30.75 * AngRad                     ' latitude of second standard parallel
Dim P0 = 29.0 * AngRad                      ' latitude of origin
Dim M0 = -84.5 * AngRad                     ' central meridian
Dim X0 = 1968500.0                          ' False easting of central meridian, map units

' Calculate the coordinate system constants.
Dim m1 = Math.Cos(P1) / Math.Sqrt(1 - (Ec2 * Math.Pow((Math.Sin(P1)), 2)))
Dim m2 = Math.Cos(P2) / Math.Sqrt(1 - (Ec2 * Math.Pow((Math.Sin(P2)), 2)))
Dim t1 = Math.Tan(Pi4 - (P1 / 2)) / Math.Pow((1 - Ec * Math.Sin(P1)) / (1 + Ec * Math.Sin(P1)), (Ec / 2))
Dim t2 = Math.Tan(Pi4 - (P2 / 2)) / Math.Pow((1 - Ec * Math.Sin(P2)) / (1 + Ec * Math.Sin(P2)), (Ec / 2))
Dim t0 = Math.Tan(Pi4 - (P0 / 2)) / Math.Pow((1 - Ec * Math.Sin(P0)) / (1 + Ec * Math.Sin(P0)), (Ec / 2))
Dim n = Math.Log(m1 / m2) / Math.Log(t1 / t2)
Dim F = m1 / (n * Math.Pow(t1, n))
Dim rho0 = A * F * Math.Pow(t0, n)
' Convert the coordinate to Latitude/Longitude.

' Calculate the Longitude.
lat = lat - X0
Dim Pi2 = Pi4 * 2
Dim rho = Math.Sqrt(Math.Pow(lat, 2) + (Math.Pow(rho0 - lng, 2)))
Dim theta = Math.Atan(lat / (rho0 - lng))
Dim t = Math.Pow(rho / (A * F), (1 / n))
Dim LonR = theta / n + M0
lat = lat + X0

' Estimate the Latitude
Dim Lat0 = Pi2 - (2 * Math.Atan(t))

' Substitute the estimate into the iterative calculation that
' converges on the correct Latitude value.
Dim part1 = (1 - (Ec * Math.Sin(Lat0))) / (1 + (Ec * Math.Sin(Lat0)))
Dim LatR = Pi2 - (2 * Math.Atan(t * Math.Pow(part1, (Ec / 2))))

Lat0 = LatR
part1 = (1 - (Ec * Math.Sin(Lat0))) / (1 + (Ec * Math.Sin(Lat0)))
LatR = Pi2 - (2 * Math.Atan(t * Math.Pow(part1, (Ec / 2))))
Debug.WriteLine(LatR & "|" & LonR)

If (Math.Abs(LatR - Lat0) > 0.000000002) Then
  ' Convert from radians to degrees.
  Dim _Lat = LatR / AngRad
  Dim _Lng = LonR / AngRad
  Return _Lng.ToString + "|" + _Lat.ToString
End If
Return ""
End Function


texts: QGIS crashes when doing a rollback after modifying values of an user-added feature with pyQGIS

texts:
I'm developing an addon for QGIS that needs to assign a UUID to any feature as soon as it is added to the layer.
I'm using the featureAdded signal to write the UUID in the corresponding field of the feature.

This is the code that I am using (as simplified as possible):

def run(self):
    self.iface.mapCanvas().currentLayer().featureAdded.connect(self.onFeatureAdded)

def onFeatureAdded(self, fid):
    layer = self.iface.mapCanvas().currentLayer()
    layer.beginEditCommand("Set UUID")
    print layer.changeAttributeValue(fid, layer.fieldNameIndex('guid_pol'), 'some_random_uuid') # prints True
    layer.endEditCommand()

(I created a new plugin with the 'Plugin Builder' add-on and that above is the only code I have added to it.)

Unless I'm missing something, this follows the guidelines from the documentation in the QGIS developer cookbook: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/vector.html#modifying-vector-layers-with-an-editing-buffer

When I edit the layer, if I add a feature the field 'guid_pol' is set to 'some_random_uuid' as expected. I can commit the changes without a problem, however, if instead of commiting the changes I try cancelling them QGIS crashes with error 0xC0000005 (Access Violation).

Interestingly, the command 'Set UUID' doesn't get added to the undo/redo stack (at least, not to the Undo/redo Panel) and if I try to undo it crashes as well.

I'm using QGIS 2.14 but I observed similar behaviour with QGIS 2.12.

I'm having a hard time getting this to work.
Is there something I am doing wrong?


texts: Is it possible to drag/draw a bbox in Openlayers 3 and have the coordinates from this function auto fill-in a form field(s) on the same page?


texts: Slope height calculation in ArcGIS Desktop?

texts:
I want to know if there exists some method to calculate the slope height from a DEM 1X1 m?

So I intend to find the rise (vertical), not the run (horizontal). I try to calculate using the slope length and the sin function of the angle, but the results are still with errors in the flat areas.

With Slope function I carried put the angle, but not the rise. So, the angle is in degrees, but the rise should be in meters, like a DEM. I didn't find any function to complete this rise calculation.


texts: How to select raster bands and layer features intersects with another layer?

texts:
I want to get raster bands from raster layer and layer features from vector layer if they intersect with another vector layer.

vector2 my main layer. If vector1 features in vector2 layer it should be selected (in this example all of them should be selected) and if raster cells in vector2 it's cell values should return.

I investigate this example but it did not work for me.


texts: ArcGIS REST API not enough storage available to complete the operation

texts:
I am building a web app for submitting jobs to the Vehicle Routing Problem Solver. I have been testing this app for a while now and I could submit jobs flawlessly before, but as of today I am getting:

ERROR 999999: Error executing function.
Not enough storage is available to complete this operation.

in the JSON response.

I am submitting the exact same job parameters every time and have done so many dozens of times before. All of the 'save' parameters for the job are set to default (false) and I have more than enough credits in my organization.

What is the cause of this issue?

EDIT: I have submitted the example job in the docs using postman and I am still getting the same error.


texts: Binding value from QComboBox in QGIS custom Form?

texts:
I am trying to build a custom form for field data collection (bird species) with some python logic. The python code inserts values into a QComboBox:

At the beginning of the list it adds the twenty most often used species
then follows a separator
then a complete list of possible species

After finding Custom QGIS feature forms – Value Binding I learned, that I have to set the correct widget type in the "Fields" tab of the layer properties dialog.

Now QLineEdits (Edit Widget -> "Text Edit") and QSpinBoxes (Edit Widget -> "Range") work as intended. But I can't find a widget type that correctly binds the values of my QComboBoxes.

Which type of "Edit Widget" do I have to choose to bind a QComboBox from a Qt custom form?


texts: Geoform for points, lines and polygons in QGIS environment

texts:
I've been asked to look for an attractive way to collect information from children and young adults about the routes they are using to go to school, sports club, etc., as well as the points on the road which they consider as dangerous.

This should be a combination between a questionnaire and a webmap on which the participants of the online survey can indicate lines (routes) and points (dangerous points).

Ideally the data of the online survey can be used afterwards in, for example, a QGIS environment for further analysis.

I know about the "Geoform", accessible through the WebApp Builder of ArcGIS Online, for which I know you certainly can collect points, but I'm looking for a free-of-charge solution for which no commercial licenses are needed.

Any suggestions about the needed configuration?


texts: Dissolving shapefile boundaries results in exogenous left over boundaries?

texts:
I'm trying to dissolve shapefile boundaries based on an attribute. I want to combine many smaller polygons so that the result is a series of polygons which are composed of the outside boundaries of the smaller ones.
I have tried to do this in R in a number of ways including:
1.
my_file <- readOGR("./my_path/my_file.shp")
groups = aggregate(my_file, by = "my_attribute")

my_file <- readOGR("./my_path/my_file.shp")
my_file_fortified <- fortify(my_file, region = "my_attribute")
# merge the "fortified" data with the data from our spatial object
my_df<- merge(my_file_fortified, my_file@data, by.x = "id", by.y = "my_attribute")

my_file <- read_sf("./my_path/my_file.shp")
st_union_by = function(geo, group) {
# browser()
geo
group
   y2 = list()
    #loop over by groups and merge units
      for (i in unique(group)) {
        #which units
            z = geo[group == i]

          #merge
          y = Reduce(st_union, z)
          y2[[i]] = y
        }

   st_sfc(y2)

}
my_file <- st_union_by(my_file$geometry, my_file$my_attribute)

taken from here: Using sf to combine polygons that share borders?
which results in this error: Error in CPL_geos_op2(op, x, y) :  Evaluation error: TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point -122.33666229315305 47.540180206469273 at -122.33666229315305 47.540180206469273.
The rest result in a map which looks like this:

Why are the broken lines bits in the polygons, and how to dissolve polygons without creating these?
This is the data and I am trying to combine based on the PROVCODE column.


texts: Unable to do a parallel INSERT using Postgres 9.6.0 & PostGIS 2.3.0

texts:
The following PostGIS query executes in parallel using Postgres 9.6's new parallel processing feature:

SELECT pnts.mb_code11, bdys.mb_code11
FROM testing.mb_random_points AS pnts
INNER JOIN testing.mb_2011_aust AS bdys ON ST_Contains(bdys.geom, pnts.geom);

Here's the query plan:

Run the same query as an insert, and it only runs on a single CPU:

INSERT INTO test
SELECT pnts.mb_code11, bdys.mb_code11
FROM testing.mb_random_points AS pnts
INNER JOIN testing.mb_2011_aust AS bdys ON ST_Contains(bdys.geom, pnts.geom);

Here's it's query plan:

How can I get the INSERT to run in parallel?


texts:
I am trying to write a VBA loop which will turn a list of addresses and their corresponding Lat Long coordinates in to a KML file for viewing in Google Earth.

I have managed to get the opening and closing KML tags in, but can't work out how to loop only the coordinates in to the coordinate tags, and have the addresses in their own tag. If anyone could lend a hand, I will be eternally in their debt! (Ultimately, I will want to have the following info in the KML File: Job Number, House Number + Street name (combined in KML Placemark tag, but from different Excel cells), City, Lat/Long. However, for my learning I currently have address and latlong.

A1: Addresses, B1: Lat, C1: Long

Current results:

"<?xml version=""1.0"" encoding=""UTF-8""?><kml xmlns=""http://earth.google.com/kml/2.2""><Document>"
"<coordinates>1892 Mccreary Road,49.950362,-97.0804002</coordinates>"
"<coordinates>38 Monarch Mews,49.948567,-97.0784119</coordinates>"
"<coordinates>3170 Vialoux Drive,49.9482365,-97.0775992</coordinates>"
"</Document></kml>"

Code:

Sub WriteTextFile()

Worksheets("Sheet1").Range("A1").Activate

Dim FilePath As String
Dim CellData As String
Dim LastCol As Long
Dim LastRow As Long
Dim CoordOpen As String
Dim CoordClose As String
Dim KMLCrap As String
Dim KMLCrapClose As String

Const HEADER As String = "<?xml version=""1.0"" encoding=""UTF-8""?><kml xmlns=""http://earth.google.com/kml/2.2"">" + ("<Document>")

LastCol = ActiveSheet.UsedRange.SpecialCells(xlCellTypeLastCell).Column
LastRow = ActiveSheet.UsedRange.SpecialCells(xlCellTypeLastCell).Row

CellData = ""

FilePath = "C:\Users\TeamCK\Desktop\VBA To KML Org\Original Examples\Great Loop.txt"

Open FilePath For Output As #2

KMLCrap = ("<?xml version=""1.0"" encoding=""UTF-8""?>") + ("<kml xmlns=""http://earth.google.com/kml/2.2"">") + ("<Document>")

Write #2, KMLCrap

CoordOpen = "<coordinates>"
CoordClose = "</coordinates>"

    For I = 1 To LastRow
        For J = 1 To LastCol
            If J = LastCol Then
                CellData = CellData + Trim(ActiveCell(I, J).Value)
            Else
                CellData = CellData + Trim(ActiveCell(I, J).Value) + ","
            End If
        Next J
            Write #2, CoordOpen + CellData + CoordClose
            CellData = ""
    Next I

KMLCrapClose = ("</Document></kml>")
Write #2, KMLCrapClose

Close #2

MsgBox ("Done")

End Sub


texts:
According to the PyQGIS Cookbook, a color ramp can be assigned in the following way:

rlayer.setColorShadingAlgorithm(QgsRasterLayer.ColorRampShader)
lst = [ QgsColorRampShader.ColorRampItem(0, QColor(0,255,0)), \
    QgsColorRampShader.ColorRampItem(255, QColor(255,255,0)) ]
fcn = rlayer.rasterShader().rasterShaderFunction()
fcn.setColorRampType(QgsColorRampShader.INTERPOLATED)
fcn.setColorRampItemList(lst)

However, this code appears to be old, because setColorShadingAlgorithm is not present in QGIS 2.2.0.  This bug report states that setDrawingStyle can be used instead, but I haven't been able to assign a color ramp through its use.

How would I assign a color ramp through PyQGIS, using setDrawingStyle or otherwise?


texts: Selenium Web Driver and OpenLayers 2.x: How to do an identify on a map?

texts:
I've to test a web mapping application that use OpenLayers 2.x, using Selenium Web Driver in Java and using Firefox (I'm on Windows 7).

I've found only this issue https://stackoverflow.com/questions/28682195/how-to-use-openlayers-drawfeature-with-selenium-webdriver-in-java-double-click that doesn't solve my problem.

I've have to test the identify function on features on the map, so:

1) select the identify button on my toolbar (I'm able to do this ... so no problem ...)

2) click on a point feature on the map (I'm not able to do this ....)

3) close the dialog that shows the felatire descriptive data (I'm not able to do this ....)

I can't give the url of my application that it's not public but I can use this simple test case

http://dev.openlayers.org/releases/OpenLayers-2.13.1/examples/markers.html

that shows my use case.

Clicking on the map, you'll see the feature details and then close the dialog.

Here you're my code that doesn't work

package myTestProjects;

import java.util.concurrent.TimeUnit;

import org.openqa.selenium.By;
import org.openqa.selenium.WebDriver;
import org.openqa.selenium.WebElement;
import org.openqa.selenium.firefox.FirefoxDriver;
import org.openqa.selenium.interactions.Actions;

public class identifyOpenLayersTest_02 {

private static WebDriver driver = null;

public static void main(String[] args)  throws InterruptedException {

    // Create a new instance of the Firefox driver
    System.out.println("Create a new instance of the Firefox driver ...");
    driver = new FirefoxDriver();

    //Put a Implicit wait, this means that any search for elements on the page could take the time the implicit wait is set for before throwing exception
    driver.manage().timeouts().implicitlyWait(10, TimeUnit.SECONDS);

    // It is always advisable to Maximize the window before performing DragNDrop action
    System.out.println("Maximize the window ...");
    driver.manage().window().maximize();
    Thread.sleep(3000L);

    // Launch the OpenLayers 2.x marker sample
    System.out.println("Launch the OpenLayers 2.x marker sample  ...");
    Thread.sleep(3000L);
    driver.get("http://dev.openlayers.org/releases/OpenLayers-2.13.1/examples/markers.html");

    // Create a new Action instance
    System.out.println("Create a new Action instance ...");
    Actions act = new Actions(driver);

    // Find the viewport inside in witch there is the map
    System.out.println("Find the viewport inside in witch there is the map ...");
    Thread.sleep(3000L);
    WebElement el = driver.findElement(By.id("OpenLayers_Map_2_OpenLayers_ViewPort"));

    // Start the action sequence
    System.out.println("Start the action sequence  ...");
    Thread.sleep(3000L);
    act.click().perform();

    // Identify marker
    System.out.println("Identify marker at 285, 111 ...");
    Thread.sleep(3000L);
    act.moveToElement(el, 285, 111).click().build().perform();

    // Print TEST = OK!!
    System.out.println("TEST = OK !!");
    //driver.quit();

        }
}

Suggestions? Samples?