texts:
I've created a Leaflet map displaying a point layer as GeoJSON. I used a custom icon (an arrow). I would like to rotate this icon on the basis of an attribute included in the GeoJSON file (field "dir"). How should I do?

JSscript code:

var map = L.map('map',{
    center: [5,28],
    zoom: 5,
    minZoom: 2,
    maxZoom: 18
});

L.tileLayer('http:...',{attribution: ...);

var icon = new L.Icon({
    iconUrl:"assets/icon/arrow.svg",
    iconSize: [30],
    });

function funzione (feature, layer) {
    layer.bindPopup("Information <br>" + feature.properties.name);
    layer.setIcon(icon);
};

L.geoJSON(data,{
    onEachFeature: funzione
}).addTo(map);

GeoJson file:

var data = {
"type": "FeatureCollection",
"crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:OGC:1.3:CRS84" } },                                                                        "features": [
{ "type": "Feature", "properties": { "time": 155011, "lat": -22.978704, "lon": -43.129650, "name": "harbour", "dir": 180} }
]};


texts: Adding an offset to a Leaflet TileLayer on the fly

texts:
I'm currently using Leaflet to display 18th century maps of Belgium as a base layer with a vector overlay of the present day road system. Although the maps are absolutely gorgeous (http://www.ngi.be/NL/NL1-4-2-3 if you're interested), map production techniques were less than ideal back then, which means some streets are seriously offset from their real positions.

What I'd like to do, is have the users drag only the 18th century base layer so they can report what offset fits a particular present day street the best.

In the following snapshot, you can see present day roads in blue, with its 18th century equivalents shifted southwards:

I'd like to temporarily drag the base layer to have the road match up somewhat to get a rough estimation of the maps (lack of) accuracy for that particular road.

It's a rather peculiar use-case, and unsurprisingly I haven't found an existing solution. I know it's possible to offset layers using Proj4Leaflet, but that seems to be at layer creation, not on the fly.

Does anyone have an idea on how to accomplish this? I'm quite well-versed in Javascript/HTML, so I'm not afraid to get my hands dirty with Leaflet internals; I currently just don't have a clue where to begin. I'd rather not move the vector layer, since there are other base layers that do match up.


texts: Submit a mass query of latitude and longitudes to soilgrid for data download

texts:
I would like to obtain soil physicochemical data for >500 sites for a study I'm doing for my dissertation research in the pacific northwest of the US. I have latitude and longitude for each site, and would like to query soilgrids with these latitude and longitudes to download the soil physical and chemical data wherever they are available.
I tried using the Explore tool on the map, and can click on individual locations and see the data are available. When I try to download a 2x2 degree tile using the download tool, I get a tif file. It would be ideal to be able to download a dataframe that includes the numerical values for all physical and chemical soil properties by latitude and longitude for multiple sites.
Is this possible with soilgrids? If not, are there any suggestions for other databases that I may query for these data instead?


texts: Shapefile string data max length and character

texts: Does anyone know why in ESRI help website (1, 2), it states that the field length if not specified will default to 255. However when I manually create a shapefile the max character is only 254. Why is there this difference of 1 character?


texts: natural=beach without all the wrong tagged sand areas golf bunker or sandy places in zoos

texts:
I would like to get natural=beach (places along a body of water at best places where you can have a bath). Unfortunately the query natural=beach ends up often with wrong "single" tagged sand areas like golf bunker e.g.: https://overpass-turbo.eu/s/sC3 or sandy places in zoos or forests. Single tagged I mean there's not any other helpful tag I could use to exclude.

Does anyone have an idea to query natural=beach without all the wrong "single" tagged sand areas like golf bunkers or sandy places in zoos or forests.


texts: Statistical method to quantify the accuracy of point clusters

texts:
I have coordinates of about 13 points from different pictures in UTM. Per point I can have 2, 3, 4 or 5 coordinates of these points from 2, 3, 4 or 5 pictures. I would like to quantify the quality of the set of pictures statistically but I don't know how. Below you find a screenshot from EXCEL to clarify the situation. If anyone knows how I could define some kind of accuracy of these images based on these point locations please let me know.

I thought about creating polygons for each cluster in QGIS (I am using 2.18) and comparing the area and width but that seems so rudimentary to me.


texts:
Is there a possible add-on that will simplify the user experience when creating a SQL statement? More specifically, is there a GUI that would allow the user to create a SQL statement similar to ArcMaps's Query Builder where pre-built function keys and data specific variables appear in the query window? CartoDB is an excellent fit for my client, however the SQL function is a little daunting to them.

I'm very new to CartoDB, but  did a fair amount of research for potential solutions without any luck.  I did researched SQL API tool, but it is my understanding that the API only interacts with the data table by Inserting, Deleting, and Updating the data?


texts: Extract pixel values using pixel condition in R

texts:
I am preparing a dataset to run a SVM classification. So far, I have a raster stack  which includes a layer ([6]) for training

> S1
class       : RasterStack
dimensions  : 3865, 6899, 26664635, 6  (nrow, ncol, ncell, nlayers)
resolution  : 14.83, 14.83  (x, y)
extent      : 361363.5, 463675.7, 5760647, 5817965  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
names       : Coherence_VV_Stack2.1, Coherence_VV_Stack2.2, Coherence_VV_Stack2.3, Coherence_VV_Stack2.4, Coherence_VV_Stack2.5, Class
min values  :                     ?,                     ?,                     ?,                     ?,                     ?,     0
max values  :                     ?,                     ?,                     ?,                     ?,                     ?,     1

Where S1[[6]] is:

> S1[[6]]
class       : RasterLayer
dimensions  : 3865, 6899, 26664635  (nrow, ncol, ncell)
resolution  : 14.83, 14.83  (x, y)
extent      : 361363.5, 463675.7, 5760647, 5817965  (xmin, xmax, ymin, ymax)
coord. ref. : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
data source : in memory
names       : Class
values      : 0, 1  (min, max)

and plot(S1[[6]) is:

I want to create a dataframe (to be used later as input for classification) where:

If pixel value in S1[[6]] is 1, extract the pixel values in the other layers of the raster stack and put it in the dataframe
If pixel value in S1[[6]] is0, do nothing.

Any suggestion on how to proceed? I know the issue can be solved by importing a shapefile and use the extract function but I want to use this approach.


texts: Arcmap: a different delimiter when symbolizing "Unique values, many"?

texts: In Arcmap 10.1, is it possible to use something other than a comma to separate the value fields when symbolizing by "Unique values, many"?


texts: Error while executing Populate_Geometry_Columns() Postgis

texts:
I loaded data into postgis using easy loader(MAPINFO). Now when I run

SELECT * FROM Populate_Geometry_Columns()

I get error :

  WARNING:  Could not convert '"SP_GEOMETRY"' in
  'vectormaps."data_bnd"' to use typmod with srid 0, type:
  MultiLineString  CONTEXT:  PL/pgSQL function
  populate_geometry_columns(boolean) line 51 at assignment

How do i check this error?


texts: Attachment Image in popup from feature layer in ArcGIS API for JavaScript 4?

texts:
I have an attachment in ArcMap attribute table and would like to show the image attachment in the popup when a feature layer element is clicked but the code I tried to show the attachment but this fails

popupTemplate.content = [{type: "attachments"}]


texts: gsb projection file for qgis under MAC OSX

texts:
I'm trying to create a new projection inside wgis. The projection includes a gsb file. With the Windows OS is easy to implemnet that: you just need to place the gsb file inside the C:\Programs\Qgis[Qgis ver]\share\proj and in Linux OS is easy as well because the gsb file goes to /usr/share/proj.

Where to put gsb files within Mac OSX?


texts: Mismatch between coordinates systems when converting from x, y -> lat long

texts:
UPDATE THIS HAS BEEN RESOLVED, THE ISSUES WERE RELATED TO THE WRONG CONVERSION OF NDSFutility JS Code into Java

I have a building's footprint which I extract from an IFC file:

https://www.iso.org/standard/51622.html.

IFC captures building geometries for buildings with a reference longitude and latitude. The building geometry is defined using a system of relative placements, and definitions of Polylines, etc. (although it supports multiple geometric representations of building geometries).

I want to convert a 2D representation of some of this geometry (walls) into geospatial data.

I have extracted the geometry (here as a list of points):

{"geometry":{"paths":[[[-2.58407192806819, -3.28072887478926], [7.715928071931811, -3.28072887478926]] ,[[-2.6840719280682, -9.080728874789271], [-2.6840719280682, -0.2807288747892702]] ,[[-13.184071928068201, -0.28072887478923203], [-13.184071928068201, -9.080728874789234]] ,[[7.7159280719318, -0.180728874789299], [-13.284071928068201, -0.180728874789299]] ,[[7.8159280719318005, -9.0807288747893], [7.8159280719318005, -0.08072887478929935]] ,[[-13.284071928068201, -9.180728874789231], [7.915928071931798, -9.180728874789231]]]}}

And I have extracted the models longitude and latitude, here represented as a wkt point POINT (-71.05977630611112 42.358428955). I want to convert the geometry to longitude and latitudes format using the given ref longitude and latitude in the IFC file so that I can display the selected 2D geometry on a map.

To do this I found this code and converted it to Java - http://www.whoi.edu/marine/ndsf/utility/NDSFutility.js

I ran the code and converted each point based on the reference longitude and latitude of the file. I converted it to WSG84. It seems that the geometries do not match though, and seem to be slightly off. Here is the wkt output:

MULTILINESTRING ((-71.05984749254343 42.358376156434254, -71.05972247208217 42.358376156434254), (-71.0599191062312 42.35832304098344, -71.0599191062312 42.35840226402872), (-71.05993974039522 42.35830773653188, -71.05993974039522 42.35822851348528), (-71.05968484441314 42.35849679152554, -71.05993974087394 42.35849679152554), (-71.05979165804781 42.3584175684816, -71.05979165804781 42.35849859204928), (-71.06004898161066 42.358226712961496, -71.05979165833584 42.358226712961496))

My code is directly mapped from the JS, so I am trying to figure out where the error is.


texts: leafR error: no applicable method for 'readLAS' applied to an object of class "c('LAS', 'Spatial')

texts:
First, I read a LAS file into R using the readLAS() function from the lidR package:
test_las <-readLAS("C:/Users/LAS files/file.las")

I am now attempting to use functions from the leafR package on the LAS file like so:
#Create a dataframe of the 3D voxels information (xyz) with Leaf Area Density values from las file
lads <- lad.voxels(file.las, grain.size = 1, k = 1)

and it produces this error:
#> Error in UseMethod("readLAS", files) :
#>   no applicable method for 'readLAS' applied to an object of class "c('LAS', 'Spatial')"

What do I need to do to prevent this error?


texts: How to apply the same style to different columns of the attributes table?

texts: How do I apply the same Categorized style, including same custom intervals  and respective labels,  to more than one column, if I have a shapefile with an extensive attributes table? If I save a style for a column, I can’t apply it to another one.


texts: How to create hillshade models every hour for a zone with GDAL

texts:
I´m trying to calculate the hillshade models from a zone for every hour during one month. I´m using the GDAL, gdaldem hillshade (I don´t know if there is another tool/library).

The problem is how to model the different situations during a day from the time and not from the elevation and azimuth of the sun. I would like use the time to generate the hillshade model for this time, so the modelling of my surface uses the different situations of elevation and azimuth for this time for my zone.

If I understood it the elevation and azimuth parameters describe the sun's position. So for the same elevation /azimuth parameters the hillshades  describe the shadows for different times (for example: the situation of the sun as "elevation 1.79 and azimuth 121" describes the shadows for a part of my zone for 01:00:00 pm and for another part for 01:30:00 pm).

Is it possible to generate hillshade models describing the shadows in a moment for a zone?


texts:
I aim to calculate the average district-specific proximity to river in India by using the proximity (raster distance) function and the following data here and the GADM India data (level 2 for districts).
First, I convert the river shapefile to a raster and use 0.0166667*0.0166667 resolution of the georeferenced units (output raster size), following these instructions. Then I apply the proximity function.
However, I am not sure what are the units of the generated distances.
If I right click on the layer it says the following:

CRS EPSG:4326 - WGS 84 - Geographic Extent
180.0000000000000000,-85.1921875000000171 : 180.0000000000000000,83.5996093750000000
Unit degrees Feature count 1,428.

Here are the proximity properties:

CRS EPSG:4326 - WGS 84 - Geographic Extent -180.0000000000000000,-85.1840615250000042 : 180.0007200000000012,83.5996093750000000
Unit degrees Width 21600 Height 10127

I tried to apply the measure tool and the values resembled mostly miles. Is this possible?


texts: Exporting 1.5 million records from File Geodatabase to Microsoft Access table?

texts:
I have a table that I want to export from arcmap that has around 1.5 million records. I need to import it into access. Now every time I have tried to export it into a format that access will read, like a csv or dbase, the records become corrupted and instead of displaying there values i they have contents such as "3.186345345345-3".

Does anyone know of a way to export to access that will work perfectly?


texts: How to set the output extents of rasterization in QGIS

texts:
I am trying to rasterize a bunch of vector based shape files in QGIS. All these files are covering certain types of ground use in a city (houses, street, open field,...).

I need these files to be rasterized because I will be running an algorithm on them to determine the least cost path. And for this I need them to be rasterized with a certain fixed pixel size (1m x 1m in the real world). My vector files are 6 km x 10 km but when I rasterize one of the files (for ex: railroads) which isn't reaching over the total 6 km in x direction and over the 10 km in y direction but only over the first 3 km x and 10 km y, the rasterized file is also only 3 km x 10 km.

Somebody any idea how to solve this? Because in this manner, when I choose a pixel size of 6000 + 10000 I'll get pixels with a size of 0,5 m x 1 m in stead of 1 m x 1 m.


texts: Is it necessary to use hibernate-spatial for postgis?

texts:
I'm beginner in GIS. I have a database in PostgreSQL that contains a table location, that table contains a column with type point. I have to get points and display it on google map. Once I did the mapping with Hibernate, it gave me for the column point the type Serializable as follow:

<property name="point" type="serializable">
        <column name="point" />

I really don't understand why it gave me serializable for that column. Is it necessary to use hibernate-spatial for that case? What can I do? I've searched a lot but I still don't know if it is necessary to use spatial hibernate knowing that I don't have much stuffs related to GIS to use for my application.


texts: Converting GeoJSON to a raster image using rasterio

texts:
My geospatial data knowledge is limited. I am hoping for advice as to how I can convert the full GeoJSON geometry into a tif image using rasterio.

I have tried using gdal_rasterize to no avail. Not all the geometry was rasterized even after I set the -at (ALL_TOUCHED) parameter.


texts: How to add wms overlay layer to google map

texts:
I have a WMS server (using MapServer and Apache on windows). i read this document about Mapserver's Tile Mode and use that's example to overlay my map on Google maps.but not work.

my code:

function load() {
  if (GBrowserIsCompatible()) {
    var urlTemplate = 'http://79.175.166.195:1280/cgi-bin/mapserv.exe?';
        urlTemplate += 'map=C:\\ms4w\\apps\\Rasht\\RashtMap.map&';
        urlTemplate += 'layers=base street&';
        urlTemplate += 'mode=tile&';
        urlTemplate += 'tilemode=gmap&';
        urlTemplate += 'tile={X}+{Y}+{Z}';
    var myLayer = new GTileLayer(null,0,18,{
                                 tileUrlTemplate:urlTemplate,
                                 isPng:true,
                                 opacity:1.0 });
    var map = new GMap2(document.getElementById("map"));
    map.addControl(new GLargeMapControl());
    map.addControl(new GMapTypeControl());
    map.setCenter(new GLatLng(37.282340, 49.588540), 15);
    map.addOverlay(new GTileLayerOverlay(myLayer));
  }
}

what is wrong?


texts: How to control visibility of cached layer in ArcGIS Server map service?

texts:
I am consuming mapservice form ArcGIS Server 10 using ArcGIS Silverlight API. I have created cache (fusion layers) for my services by right clicking on the mapservice in ArcCatalog and clicking create cache. It works fine in my application but the problem is after viewing the cache when do map full extent the cache is still visible when it had to hide since I have set zoom visibility.

Please guide me as to how I can keep the map cache visibility off when it is out of the zoom visibility.


texts: How is browsers position altitude defined?

texts:
Looking at the docs:

https://developer.mozilla.org/en-US/docs/Web/API/Coordinates/altitude

It says "altitude above sea level in meters". Does this mean it takes into account geoid data of some kind?


texts: eCognition TTA Mask : Aerial images in ecognition for supervised classification

texts:
I am using eCognition Developer Trial 8.0. I have rgb Drone Otho foto. We want to classify this image for various biotope types exist in that region. We have already an existing classification as a thematic layer that contains at least 100 of classes.

My primary idea was to classify based on TTA mask. So that I donot need to assign it by hand. I have genererated TTA mask using this solution.

I do not have any idea if it is possible to classify the whole area/image based on TTA mask or how it works.
You can take a look of my data screenshots:


texts: How do I convert affine coordinates to lat/lng?

texts:
I am extremely new to GIS.

I'm using gdal to read in a landuse/landcover map and I need to pick out the lat/lng of certain land cover types to index into a different dataset which is expressed only in lat/lng.  Unfortuantely, I don't understand the form of the x and y coordinates given to me from the geotransform, specifically the originX and originY below:

geotransform = dataset.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]

Printing these values gives me coordinates like (447466.693808, 4952570.40529).  How do these relate to the original latitude and longitude?

Edit:

Here's a simple python example that got me what I was looking for:

srs = osr.SpatialReference()
srs.ImportFromWkt(dataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs,srsLatLong)
print ct.TransformPoint(originX,originY)

Stolen from: tolatlong.py


texts: Integrate ArcGIS 10 extension INTO ArcEngine 10

texts:
I was wondering if there is a way at all to use some of ArcGIS 10's extension into arc engine.

A good example would be if we could put Arc Hydro extension into ArcEngine.


texts: How do i identify unknown coordinates reference system for my data

texts: I have  several shapesfiles which perfectly overlay but the coordinates reference system is unknown. I need to overly these with other data in a different coordinate system.  I tried initially to overlay with several projected and geographic coordinates systems associated with the area out of which two showed very close displacement (about 60m). Next i tried different datum with one of the 'promising' coordinates reference systems (a geographic coordinates system) and also changed the units but same results. Attached is a clip of one of the coordinate extents. I am yet to try the second 'promising' CRS but decided to put it out just in case i am missing something


texts: In Qgis, how do I move a group of point/attributes with an offset of 100 metres

texts:
I have been using Qgis for a few years now, not high level stuff, and not professionally trained.

I have a set of attributes on a survey from a plane, and would like to move the entire group of points, lets say 100 metres East or West from where the animals are. By being able to offset the right side of the plane count and the left side of the plane count it is easier to look at and make sense of. The survey is all about numbers not their specific location.

The screen shot is a ground based test but shows green symbols being the right side of a plane or vehicle and the purple the left side (there were two runs on the left)

How can I do this?


texts:
I want to draw multiple buffer selection on click so that each selection is saved in a different vector layer and buffer features in one vector layer. The following code works well for the first buffer and selection layer. But on second time two empty selection are created and 3rd time 3 selection empty layers and so on, what I could figure out is that featureadded event for the buffer layer is called no of times of the total features in bufferLayer every time.

map.addLayers([bufferLayer, pointLayer]);
pointCtrl = new OpenLayers.Control.DrawFeature(
        pointLayer, OpenLayers.Handler.Point, {
        eventListeners : {
            //beforefeatureadded    : function (e) {e.attributes = { id: featureId++, name : "Buffer Point "+featureId, radius: radius, sides:sides }},
            featureadded        : function (e) {
                var center_point = new OpenLayers.Geometry.Point(e.feature.geometry.x, e.feature.geometry.y);
                var circle = new OpenLayers.Geometry.Polygon.createRegularPolygon(center_point, radius, sides, 0);
                var circleFeature = new OpenLayers.Feature.Vector(circle);
                circleFeature.attributes = { id: featureId, name : "Buffer "+featureId, radius: radius, sides:sides };
                bufferLayer.events.on({
                    featureadded        : function(e){
                        var x;
                        var geometry = e.feature.geometry.clone();
                        var spatialfilter = new OpenLayers.Filter.Spatial({type : OpenLayers.Filter.Spatial.INTERSECTS,value : geometry.transform(mercator, geographic)});
                        bufferFilter =      new OpenLayers.Filter.Logical({type : OpenLayers.Filter.Logical.AND, filters : [spatialfilter, getGeneralFilter()]});
                        var x;
                        var cloneLayer=hfLayer.clone();
                        cloneLayer.filter = bufferFilter;
                        cloneLayer.name="Buffer Selection "+featureId++;
                        console.log(featureId);
                        cloneLayer.refresh({force : true});
                        cloneLayer.events.on(
                                {
                                featuresadded : function (event) {
                                    pointCtrl.deactivate();
                                    pointLayer.styleMap.styles.default.rules[0].symbolizer.label = "Health Facilities\n Within " + radius / 1000 + " Kilometers\n" + cloneLayer.features.length;
                                    map.zoomToExtent(bufferLayer.getDataExtent());
                                }, scope: this
                                });
                        cloneLayer.refresh({force : true});
                        map.addLayers([cloneLayer]);
                    },
                        scope : this
                });
                bufferLayer.addFeatures(circleFeature);
                pointCtrl.deactivate();

            },
            scope : this
        }
    });
map.addControl(pointCtrl);
pointCtrl.activate();


texts: Is it possible to get OSM data in geojson via overpass api using ipython?

texts:
I am looking for a solution that works with the overpass api.

for example:

  http://overpass-api.de/api/interpreter?data=[out:geojson];area(3600225494);(node(area)[%22amenity%22~%22theatre|cinema|arts_centre%22];);out%20body;

When trying this approach I get an error that says

  "Error: line 1: static error: For the attribute "output" of the
  element "osm-script" the only allowed values are "xml", "json", "csv",
  "custom", or "popup"."

Therefore GeoJSON is not supported this way.

I tried several approaches in ipython (python 2.7) to transform json into geojson, but failed so far.

What I need is a way to get OSM as geojson in ipython via using a url (e.g. overpass api).

BR


texts: How do I use actions from ArcGIS SDK for WPF in the latest version of WPF?

texts:
I understand how to use the Actions and Behaviors but there is a problem I'm not sure is solvable other than to downgrade the version of .net and WPF I'm using.

I installed the latest ArcGIS SDK for WPF, version 10.2 and saw one of their examples illustrating how to use the MeasureAction. When I try to implement this in my own WPF program it won't compile and tells me that the version of Blend that MeasureAction was compiled with is different than the one I'm using. I am using the most recent Microsoft.Expression.Interactions and System.Windows.Interactivity (file version says 4.0; NOT assembly version) but in their examples they are using version 2.0 (again, file version not assembly version).

Why would Esri release a new SDK (fairly recently) and not compile it with the most current version of those assemblies? Is there a way to make this work in my project? I can't roll back to a different version of .net because I'm using Caliburn.Micro which requires the version I'm using now.


texts: QGIS GRASS v.distance function not producing results

texts:
I am trying to use the v.distance function from QGIS ToolBox to calculate the closest distance between points to areas. The problem is that even if I have added new columns ("dist", "attribute") to the from_point file (.csv), there are no options for me to choose them for the upload column (using the GRASS tool GUI; please see the screenshot below). The only options are the lat long columns, which obviously means my point coordinates get overwritten immediately and the processing runs into an infinite loop of errors...
I am now trying to look into running the grass function on Python Console, but it hasn't worked so far (no error message but no output either; the upload column is empty after running).


texts:
I have a very large shapefile of building polygons and a shapefile of the country they reside in (England). The building shapefile was built by merging 94 smaller shapefiles. Both files have the same CRS.

I have been trying to use the difference tool in QGIS to create a country polygon with holes in where the buildings are.

Oddly the tool completes after about an hour and adds the new difference layer but it is blank. I have 12 GB of RAM and 64bit QGIS but still suspect it's a memory/software issue as I tried the exact same worfklow with just a 25 sq km section of the polygons and it was fine.

Can anyone recommend alternative software (or alternative workflows) to try?

Examples: I wondered if I could load it into a PostGIS database and do a spatial query for example? OR if i could create a batch process to calculate the difference for each building file one after another.


texts: Creating a buffer based on the number of parallel lines

texts:
I am working on a project for a gas pipeline and I am trying to create a study area based on the pipeline alignment. The problem is that there are multiple pipelines within the right-of-way and the study area must be different widths based on the number of lines that run parallel to one another. For example, a single pipeline must have a 25' buffer, a double pipeline must have a 37.5' buffer and a triple pipeline must have a 50' buffer.

Is it possible to do this without going through and manually selecting the linework and individually buffering them?

I have been working on this for a few hours and I have had no luck.

I am using ArcMap and I have an advanced license.


texts: Applying geometric operation on multiple GeoDataFrames

texts:
I have a collection of geodataframes in my Python environment. How can I apply for each GeoDataFrame the same spatial operation using following manner and create a new variable with the reference that it was manipulated by attaching "_dissolved"?
shp_1 = gpd.read_file('test1.shp')
shp_2 = gpd.read_file('test2.shp')
shp_3 = gpd.read_file('test3.shp')

lst = [shp_1,shp_2,shp_3]

for i in lst:

    i = i.dissolve(...)

Expected result would be the variables:
shp_1_dissolved
shp_2_dissolved
shp_3_dissolved

I tried something like:
i'_dissolved' = i.dissolve(...)

or
i + '_dissolved' = i.dissolve(...)


texts: Loading spatialite in Python (mod_spatialite vs. libspatialite.so)

texts:
I would like to write Python scripts that can load Spatialite across Linux platforms and Spatialite versions. For Spatialite versions lower than 4.2, I can use this:

from pysqlite2 import dbapi2 as sql
con = sql.connect("mydb.sqlite")
con.enable_load_extension(True)
con.execute('SELECT load_extension("libspatialite.so");')

For Spatialite versions >= 4.2, I must change the last line to:

con.execute('SELECT load_extension("mod_spatialite");')

Neither is cross compatible for me. Is there a solution for cross compatibility I am missing? I could invoke try and except using both load_extension commands, but I would prefer a one-line solution. Does one exist?


texts: Calculating TOA reflectance from Radiance using Erdas Imagine for a Spot2 image

texts:
Using the SPOT equation   ρ = π L/[E cos (θs) d-²] to calculate TOA reflectance from Radiance.

L was calculated using this equation [L = (X/A) + B]

d was calculated by converting date to julian day then converting to d using http://landsathandbook.gsfc.nasa.gov/pdfs/Landsat_Calibration_Summary_RSE.pdf

Esun  was calculated using http://www2.astrium-geo.com/files/pmedia/public/r452_9_normalsolarirradiance.pdf

To calculate cos (θs), taking the sun elevation from 90 and then converting to radians http://www.mathinary.com/degrees_radians.jsp  and then getting the cos of the answer

Example of above equation filled in for a Spot2 image p= (3.14*radian band1)/(1865*(0.98717*0.98717)*cos(70))

This equation works in Erdas Model maker however the resulting image is blank.

I need to know where I'm going wrong using the equation in Erdas


texts: Function process diagram of cut polygon tool

texts:
As far as the Cut Polygon tool from the Editor toolbar is concerned, polygons can be cut into two or more shapes according to that line drawn by a user.

Where can I get a detailed function process diagram of that tool, or any other tool or command of ArcGIS Desktop?


texts:
Code:

 Dim name_array() As String
 Redim name_array(10)

How to add value to the array? I know we can add value by specifying index like   name_array(i) = "some_value"

But I don't want to specify any index i. Instead is it possible to Add a value to array to last index used.

eg.  If i have added to i=5 Will I be able to add to 6th index without specifying index?


texts:
I'm trying to reduce the amount of physical field pickup for a feature survey plan by flying the property, tracing over the point-cloud and exporting the linework out as a DXF so I can edit it in my CAD software (12D model, not a huge following for it's drafting ability, very good for survey work though)

I've been able to do the following;

Create a custom CRS for the map projection that we use for the bulk of our survey work and what Photoscan also uses
Import the georeferenced orthometric photo as a raster layer
Import the georeferenced DEM as a raster layer and set transparency to 50% so I can see my "base map" underneath
Create and draw lines by Shapefile Layer (even if I don't fully get the fields list at this stage)
Export the Shapefile as a DXF, with correct CRS (I know it's right as when I view it on my CAD software it sits exactly as it should, 2D anyway) with Geometry type with several types (LineString/Polygon/GeometryCollection) with "Include z-dimension" ticked

I can import it in fine, 2D sits fine, but there is no height attached. I am, probably naively, thinking that since I'm tracing over the DEM then there should be z-values attached but obviously there isn't.

Is there something I'm missing out on?


texts: Classifying areas in scanned maps (not by colors but by pattern)

texts:
I have old land-use maps in a scanned format. I am trying to obtain proper land-use maps from it, similar to what is decribed here: Comparing areas in scanned historical land use maps using QGIS?.
The problem I'm facing is that some of the land-use classes are drawn as pattern (e.g. forests below are green with sprinkled black dots). For many techniques in GIMP, such as selecting areas by colors / selecting all green pixels, which are in close to black pixels, the issue is the dithering and the labels/raster of the map. The colors are fairly easily recognized with any software.

My fairly general question is for tips on how to do this in ArcGIS/GiMP? I.e. possibly with the intermediate step of creating an image file, where the patterns are replaced by plain areas in high-contrast colors.

So far I've tried colour-filters and various ways of "selecting by color" in gimp.

Example from the map:


texts: How to calculate multiple variables with precondition

texts:
I want to do multiple "calculate values" in ModelBuilder.
To make it simple let´s take only 2 simple calculations:
1) calculate sine from angle to variable "MySinus"
2) use this value in a second calculation to calculate a hypothenuse in a triangle.

I want to keep the result of calculation 1) in a separate Variable "MySinus", so doing 1) and 2) in a single step is not an option (the real calculation i have to do is more complicated with 15 and more connected calculations plus geoprocessings)

It is clear, that calculation 1) must be done before calculation 2).
I could use a precondition, but what if the result is a negative number? Then the 2nd calculation won´t be executed, because precondition is "false".

Calculation Value 1:
Expresssion: math.sin(-math.pi*0.55)

Calculate Value 2:
Expression: %abc%*5

What I tried already:
- Calculation Value 1) data type to "variant"
with result:
Executing (Calculate Value 1): CalculateValue math.sin(-math.pi*0.55) # "Any value"
Start Time: Mon Jan 12 11:39:10 2015
Value = -0,987688340595138
Succeeded at Mon Jan 12 11:39:10 2015 (Elapsed Time: 0,28 seconds)
Executing (Calculate Value 2): CalculateValue -0,987688340595138*5 # Double
Start Time: Mon Jan 12 11:39:10 2015
Value =
Succeeded at Mon Jan 12 11:39:10 2015 (Elapsed Time: 0,05 seconds)

(Second Value isn´t really calculated !!!)

Calculation Value 1) data type double
Result:
Executing (Calculate Value 1): CalculateValue math.sin(-math.pi*0.55) # Double
Start Time: Mon Jan 12 11:40:50 2015
Value = -0,987688340595138
Succeeded at Mon Jan 12 11:40:50 2015 (Elapsed Time: 0,05 seconds)
Executing (Calculate Value 2): CalculateValue -0.987688340595138*5 # Double
Start Time: Mon Jan 12 11:40:50 2015
The process did not execute because the precondition is false.
Succeeded at Mon Jan 12 11:40:50 2015 (Elapsed Time: 0,00 seconds)
Change Expression in Calculation 2) to: float(%abc%)*5
Result:
Executing (Calculate Value 2): CalculateValue float(-0,987688340595138)*5 # Double
Start Time: Mon Jan 12 11:42:12 2015
ERROR 000539: Error running expression: float(-0,987688340595138)*5
Traceback (most recent call last):
File "", line 1, in
TypeError: float() takes at most 1 argument (2 given)

Failed to execute (Calculate Value 2).
Failed at Mon Jan 12 11:42:12 2015 (Elapsed Time: 0,05 seconds)

(it seems that there is a problem with comma recognized as parameter delimiter)

Change Expression in Calculation 2) to: float("%abc%".replace(",","."))*5

Executing (Calculate Value 1): CalculateValue math.sin(-math.pi*0.55) # Variant
Start Time: Mon Jan 12 11:53:47 2015
Value = -0,987688340595138
Succeeded at Mon Jan 12 11:53:47 2015 (Elapsed Time: 0,07 seconds)
Executing (Calculate Value 2): CalculateValue float("-0,987688340595138".replace(",","."))*5 # Double
Start Time: Mon Jan 12 11:53:47 2015
Value = -4,93844170297569
Succeeded at Mon Jan 12 11:53:47 2015 (Elapsed Time: 0,05 seconds)

Finally !!!!!

Is there any simplier way instead of this there an back conversion??
(Of course I could use a python script for all this calculations, which seams to be straight way. But how to do it easy in ModelBuilder?)


texts: How can I do a subquery inside a subquery in PostGIS?


texts:
I am using OpenLayers API in my web application and trying to interact with the local GeoServer via WPS. The problem I am facing is that, WPS client.execute() works only in case, if the output of the Process is Geometry type (i.e. forJTS:intersection,JTS:buffer, etc). It does not seem to work if the output is double/string/bool(i.e. forJTS:area,JTS:envelope,JTS:equalsExact, etc).

For e.g., the following code does not popup a "success" message.

 var geometry = OpenLayers.Geometry.fromWKT(
                   'POLYGON((110 20,120 20,120 10,110 10,110 20))'
                );

 client.execute({
            server: "local",
            //process: "JTS:envelope",
            process: "JTS:geometryType",
            inputs: {
                geom: geometry
            },
            success: function(output) {
                alert('success');
            }
  });

But if I happen to replace JTS:geometryType with JTS:envelope, it works !!..

Not much of tutorials are available on internet. Is it a restriction in client.execute() or am I doing something wrong in the code ?


texts: Carto DB not using all available CPU resources

texts:
I have installed a Carto instance using the Carto docker container found here using the installation process outlined at that link (nothing more, nothing less). When performing normal functions, only a single thread is utilized causing large geometries and/or analysis to take a while to load. Is there any way to increase utilization, or am I stuck with a single threaded process?
Here's an example:

And it only uses one thread:

It will run for a while (15+ mins), then finally load... When working with large datasets (1+ gb worth of kml files), it makes it difficult to impossible to create usable maps...


texts:
I was wondering if there was a quick way to merge multiple .tab files together? I have about 100.

Would appending one tab to another with the MapBasic window open, then changing the name of the table appended to TableInfo(1,1) work?


texts:
I have two shapefiles, one polygons the other points.

I'm using Data Driven Page Setup, to produce a series of layout.

In the Index Layer I chose the shapefile of the polygons that have the OBJECTID that I use as reference for each Polygon.

How can I add the shapefile of the points so that I can get the X, Y, N_Order information that is in the Attributes table to get this information in the Dynamic Layout?


texts:
I want to plot a raster of the forest cover accross Europe using tmap::tm_raster(). The raster data are available here:
If I plot it, this is what I get:

# Read raster data:
forest <- raster("eea_r_3035_100_m_forest-area-2015_p_2015_v1_r1.tif",      NAvalue = 65535)

# Plot data:
  tm_shape(forest_br) +
    tm_raster(col = "eea_r_3035_100_m_forest-area-2015_p_2015_v1_r1.tif",
              palette = "darkgreen",
              title = "Forest cover")

when I check the raster properties, I see that I have an NA value equal 65535:
> forest
class      : RasterLayer
dimensions : 46000, 65000, 2.99e+09  (nrow, ncol, ncell)
resolution : 100, 100  (x, y)
extent     : 9e+05, 7400000, 9e+05, 5500000  (xmin, xmax, ymin, ymax)
crs        : +proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
source     : C:/input/eea_r_3035_100_m_forest-area-2015_p_2015_v1_r1.tif
names      : eea_r_3035_100_m_forest.area.2015_p_2015_v1_r1
values     : 0, 65535  (min, max)

My forest = 1, land = 0, and NA should be transparent. If I set NAvalues(forest)<-65535, my min-max values remains the same.
Please, how can I specify this in tmap??? palette seems to works only for continuous data, but what about in my raster is binary and have many NA values (resp. 65535)?
THis is how it should look like (raster reclassified outside of R. Reclassify the values using classical matrix takes forever on my Win 10 64 b, R version 4.0.2 (2020-06-22)).


texts: Identifying a layer in qgis from clicking the canvas

texts:
I have a problem with identifying shapefiles in the layer list. I am working on one project with about 50 different shapefiles representing fishing areas placed upon a nautical map. I don’t know what all of them are called by head so I would like to be able to click on them in the map canvas and somehow see which one it is in the layers section.

I hope I could explain the problem clearly. I have tried to find a tool for this but not successful. There must be a simple solution I hope.


texts: How can I display popup when graphic object is added in map using Arcgis Javascript

texts:
We have to display callout popup for each graphic object that are available on map. i.e if point is created on map, it should be associated with callout popup which will display information about that point graphic. No need to create on point to show popup. So if 5 point object are on Map, callout for all 5 point should be displayed at once when these are created.

My question is how can I display such popup always for all points that are created on map?


texts: Clipping/cropping/masking a raster tile layer in Mapbox GL JS

texts:
I'm loading Landsat tiles into a Mapbox GL JS map via an AWS S3 tile server, and this is working OK, but what I'd really like to do is clip the tile layer to a GeoJSON polygon which is loaded in another layer. By this I mean instead of showing the whole image, which in this case covers a huge area, it would just show the part of the image that lies within the boundaries of the GeoJSON polygon - in this case an area of approx 2.5ha. Note that I'm working with irregular polygons too rather than just rectangular bounding boxes.

My first thought was that I might be able to achieve this by somehow setting the raster layer in a fill parameter for the GeoJSON layer, but that doesn't seem to be possible. Then I thought maybe the GeoJSON layer could be given as a 'mask' or 'clip' parameter when loading the raster layer, but there doesn't seem to be any support for this either.

Has anyone else worked out a way to do this? I'd really like to do it with Javascript if possible, but Python may be an option.


texts: Difficulty drawing map lines across the international date line; Connection Map Visualization in R

texts:
I'm creating a visualization from country centroids to a few thousand coordinate pairs, and I'd like them to all be displayed unbroken and smooth on the same map, but I get these horizontal lines when passing the international date line at the edges. I'm using a modified version of this nycflights13 visual. My error is the same as this GIS Stack Exchange user's error, but using R.

My code is below:

library(maps)
library(geosphere)

df <- data.frame(
  Lat_pt = c(24.30515,23.86345,23.39825,24.17565,23.09845,21.46860),
  Lon_pt = c(85.14140,89.73296,81.40266,88.11363,85.85935,90.96883),
  Lat_centroid = 45.67626,
  Lon_centroid = -112.4708)

maps::map("world", fill = T, col = "grey8", bg = "grey15")
points(df$Lon_pt,df$Lat_pt, pch=3, cex=0.1, col="chocolate1")

for (i in (1:dim(df)[1])) {
  inter <- gcIntermediate(c(df$Lon_pt[i], df$Lat_pt[i]), c(df$Lon_centroid[i], df$Lat_centroid[i]), n=200)
  lines(inter, lwd = 0.5, col="turquoise2")
}

Which generates this visual:

How can I force these paths that draw across the date line to instead draw only on the map, directly from centroid to point?

I believe my error results from my use of the geosphere package, documented here


texts: Reading in AW3D30 GeoTIFF files and processing them in Python?

texts:
I am new to GIS processing, and I am using Python 3.6 to do so.
My current task is to read in several files of geographic data (specifically tiles of elevation data from JAXA), combine them for my area of interest, and then analyze them using Python.

Unfortunately, the data for each tile from AW3D30 come as a folder of six files: three tif files and two txt files.  I haven't found any documentation about what each of these files contains, so I don't know which one(s) I need for the elevation data (although it must be at least the DSM file).  And I don't know how to properly read them in and process them together, connect the proper lon/lat data to the pixels, etc.

For example, I read the *_AVE_DSM.tif file into Python using gdal, found out there is one band, and the band has a min and max value of None.

altitudeData1 = gdal.Open('N034E138/N034E138_AVE_DSM.tif')
print("[ RASTER BAND COUNT ]: ", altitudeData1.RasterCount)       #-> 1
print("[ MIN ] = ", altitudeData1.GetRasterBand(1).GetMinimum())  #-> None
print("[ MAX ] = ", altitudeData1.GetRasterBand(1).GetMaximum())  #-> None

So I'm clearly missing something important.

My final desired output is geoPandas dataframe with a geometry column capturing the polygon of each pixel (i.e., not the center point, but the square's geometry) and a column for elevation.  I am open to various workflows and intermediary steps.


texts: Pansharpening Sentinel-2 imagery in Google Earth Engine

texts:
I am attempting to pansharpen Sentinel-2 imagery in GEE in order to determine the amount of land covered in water at my field site locations (probably will calculate MNDWI on sharpened image). I do not have any experience in GEE or with pansharpening, so I am not sure whether my issue is a coding issue or a mathematical issue. I have been following this code < https://leclab.wixsite.com/spatial/post/pansharpening-sentinel-2-imagery-in-google-earth-engine >, which pansharpens the Sentinel imagery to all 10m bands using a high-pass filter (HPF) resolution merge. I have been following the code using my own Sentinel-2 scene, but am receiving an "internal error" when I try to print the result of the MSstdev calculation in Section 5c. I can print all images and numbers before this step, and nothing after this step will print (although the entire code will run without any errors).

I have included a snippet of the code below where the error occurs, if anyone might know what the issue is.


texts: Shapefile postgis database to be queried by GeoDjango

texts:
I would like to create a Postgres database (postgreslq-9.1 and postgis-2.0.2) storing several shapefiles but I am uncertain on how to implement it so I'm searching for some advices.

I have already done this last year but was pretty uncofortable beacause I had a table per shapefile, every table signed with its SRID (stored inside geometry_columns table - I was using postgis 1.5, so this was still a table and not a view as in postgis-2).

I had to create a new table every time I converted a shapefile: pretty uncofortable because I had to query them using Django ORM, and create table mapping on the fly is not trivial.

Is that possible to create one single table having an id indicating the original shapefile and the relative SRID per geometry?

Otherwise have you got any better ideas?

Thank you.


texts: Calculating areas of different raster classes in R with Superclass

texts:
I have a classified image with three categories (specific grass species: 1, ground: 2, other vegetation: 3). I am trying to calculate the area of each category. For the classification process I used the superClass RStoolbox method. For the calculating area method, I am using another code in a similar question, but am having some trouble with my image being in a superClass format type. I tried exporting and importing the file as a new raster and ended up getting a different error. Any suggestions?

img <- brick("C:/Users/name/location/Blung/null.tif")
shp <- st_read("C:/Users/name/location/Blung/Training.shp")
shpS4 <- as_Spatial(shp)

NDVI.Overlay <- function(b1, b4) {
  NDVI.Calc <- (b4 - b1) / (b4 + b4)
  return(NDVI.Calc)
}

NDVI <- overlay(img[[1]], img[[4]], fun = NDVI.Overlay)

img_update <- addLayer(img, NDVI)
names(img_update) <- c('b1', 'b2', 'b3', 'b4', 'NDVI')

SC <- superClass(img, trainData = shpS4, responseCol = "id", model = "rf", tuneLength = 1, trainPartition = 0.8)

SC[] = sample(1:3, ncell(SC), replace=TRUE)
counts <- tapply(area(SC), SC[], sum)

Error in (function (classes, fdef, mtable)  :
  unable to find an inherited method for function ‘area’ for signature ‘"superClass"’

writeRaster(SC$map, filename = "C:/Users/name/location/Blung_classification_R.tif")
img_processing <- raster("C:/Users/name/location/Blung/Blung_classification_R.tif")

img_processing[] = sample(1:3, ncell(img_processing), replace=TRUE)
counts <- tapply(area(img_processing), img_processing[], sum)

Error in res[i] <- readBin(x@file@con, what = dtype, n = 1, size = dsize,  :
  replacement has length zero
In addition: Warning messages:
1: In .local(x, ...) :
  This function is only useful for Raster* objects with a longitude/latitude coordinates
2: In .rasterFromRasterFile(grdfile, band = band, objecttype, ...) :
  size of values file does not match the number of cells (given the data type)


texts: Managing intermittent data in ModelBuilder?

texts:
I'm writing a model to automate combining county land cover masks in a large watershed. For example, I run Mosaic to New Raster to combine a Rail Right of Way mask (1,null) with an impervious structure mask, and an impervious Non-Road mask. Many Counties have Rail Land Cover and many others do not. Even if the Rail does not exist, I still want to mosaic any specified masks that do exist.

I know I can use "Exists" and Merge Branch, but I actually have many different classes that do not always appear in a given county, e.g. beaches, tidal wetlands, and others. It becomes onerous to write Calculate Value and Merge Branch conditions for every permutation of a set of classes existing or not. I see no way to create a "Null" or NoData Raster as dummy input.

Creating Zero rasters for each missing datatype involves huge amounts of disk space as I am working at 1 meter resolution in a very large watershed. Any idea how I can get Mosaic to Raster to go ahead and Mosaic all the existing raster input I specify and ignore those that do not exist? Does anyone know of a more or less elegant work around?


texts: Set projection (netCDF) with R-ncdf4-Package

texts:
I want to create a netCDF-file that has a defined projection ("WGS 84 / UTM zone 32N") with the netcdf4 package in R. The input float array has the size of x=880 and y=760. So i define the two dimensions:

ext = extent(raster_stack.ras)
res = 50
x = seq(ext[1]+res/2,ext[2]-res/2,50)
y = seq(ext[3]+res/2,ext[4]-res/2,50)

dimx = ncdim_def(name = 'x', longname = 'x coordinate of projection', units='m', vals = as.double(x))
dimy = ncdim_def(name = 'y', longname = 'y coordinate of projection', units='m', vals = as.double(y))

The extent gives the xmin, xmax, ymin and ymax of the input raster and define the columns and rows of the array. To get the midpoint of one pixel it's necessary to add a half of the pixelsize (final pixelsize should 50 m). Finally with dimx and dimy the dimensions of the netCDF are getting produced.

Then I define the variable as follows:

varz = ncvar_def(name = out.layername, longname = out.layername, units = 'percent', dim = list(dimx,dimy), missval = -99, prec = 'float')

outnc = nc_create(out.name,varz,force_v4=TRUE)
ncvar_put(outnc,varz,soil_moisture)
nc_close(outnc)

If I create the netCDF there is no projection for the file. But i can't find any documentation how to define the projection or set it in a variable.

To have an example for the netCDF Header, I converted a geoTIFF with gdal_translate to netCDF and get the following header. There you could see the defined projection. How could I declare it in R with netCDF4 ?

 2 variables (excluding dimension variables):
    char transverse_mercator[]
        grid_mapping_name: transverse_mercator
        longitude_of_central_meridian: 9
        false_easting: 5e+05
        false_northing: 0
        latitude_of_projection_origin: 0
        scale_factor_at_central_meridian: 0.9996
        longitude_of_prime_meridian: 0
        semi_major_axis: 6378137
        inverse_flattening: 298.257223563
        spatial_ref: PROJCS["WGS 84 / UTM zone 32N",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["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",9],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["metre",1,AUTHORITY["EPSG","9001"]],AUTHORITY["EPSG","32632"]]
        GeoTransform: 550000 50 0 5740000 0 -50
    float Band1[x,y]
        long_name: GDAL Band Number 1
        _FillValue: -99
        grid_mapping: transverse_mercator
 2 dimensions:
    x  Size:880
        standard_name: projection_x_coordinate
        long_name: x coordinate of projection
        units: m
    y  Size:760
        standard_name: projection_y_coordinate
        long_name: y coordinate of projection
        units: m


texts: How can I reclassify a raster to floating point values?

texts:
I have a raster of landuse values that I need to reclassify into pollution concentration values. I have the pollution values in a separate table. They are decimal values (ex. 0.0056), and my landuse values are integers 1 through 6. I have a table mapping the landuse values to the pollution values.

I tried using the Reclass By Table tool, setting my pollution values as the "Output value field", but the tool description says that this field must be an integer field. How can I reclassify the values into decimal values? I feel like I must be missing a much easier, alternate solution.

My goal is to have a pollution concentration raster that I can use later in further calculations.


texts:
When trying to use the ZoomToPostcode plugin in QGIS 2.16.2 Nødebo, after inputting a postcode to jump to I am prompted to download updated postcode files. Clicking Yes generates a "Unable to download file" dialog, and nothing else.

The Python Console outputs the following (username redacted for privacy):

File "C:/Users/xxx/.qgis2/python/plugins\ZoomToPostcode\zoomtopostcode.py", line 138, in check_pkl
    self.update_pcode_option()

File "C:/Users/xxx/.qgis2/python/plugins\ZoomToPostcode\zoomtopostcode.py", line 176, in update_pcode_option
    self.download_pkl()

File "C:/Users/xxx/.qgis2/python/plugins\ZoomToPostcode\zoomtopostcode.py", line 209, in download_pkl
    pcode_zip = zipfile.ZipFile(zip_temp)

UnboundLocalError: local variable 'zip_temp' referenced before assignment

How can I resolve this, and have a functioning postcode search function?


texts:
I use QGIS for GIS. Lately I had just too many points appearing on my vector layer and I will want to reduce them (100 for a given scale : min-1 and max 9000). With labels it works quite well but how can this be done with the features on the vector layer in QGIS?
I will need an expression to choose at random 100 points for the layer, that is if the query builder could do this.


texts: How to render layer into png in a standalone application?

texts:
I would like to use standalone appliction to apply .qml style to a layer and use QgsMapRenderer() to save it as image on a hard disk. The are a lot of questions like mine at this site and it seems that most of the issues emerge from settings for QgsApplication(). Unfortunately it seems that I wasn't able to setup QgsApplication() in a right way. My code is pretty standart rendering code from PyQGIS Cookbook that was used as a template in other questions like mine. It produces transparent .png inage with no layer rendered. Though this code works (without initialisation of QgsApplication()) if I ran it line by line in QGIS console (QGIS 2.4).

So I have following QgsApplication() settings (on Alt Linux).

Application state:
QGIS_PREFIX_PATH env var:
Prefix:     /usr
Plugin Path:        /usr/lib64/qgis
Package Data Path:  /usr/share/qgis
Active Theme Name:
Active Theme Path:  :/images/themes//
Default Theme Path: :/images/themes/default/
SVG Search Paths:   /usr/share/qgis/svg/
User DB Path:   /usr/share/qgis/resources/qgis.db

The code itself:

from qgis.core import *
from qgis.utils import *
from qgis.gui import *
from PyQt4.QtCore import QSize
from PyQt4.QtGui import QPainter, QImage, QColor

QgsApplication.setPrefixPath('/usr', True)
QgsApplication.initQgis()
print QgsApplication.showSettings()

lyr = QgsVectorLayer('test.shp', 'points', 'ogr')
lyr_id = lyr.id()
QgsMapLayerRegistry.instance().addMapLayer(lyr)
lyr.loadNamedStyle( 'test.qml')

# create image
img = QImage(QSize(800, 600), QImage.Format_ARGB32_Premultiplied)

# set image's background color
# color = QColor(255, 255, 255)
# img.fill(color.rgb())

# create painter
p = QPainter()
p.begin(img)
p.setRenderHint(QPainter.Antialiasing)

layers = QgsMapLayerRegistry.instance().mapLayers()
lst = layers.keys()
print 'lst', lst
renderer = QgsMapRenderer()
renderer.setLayerSet(lst)

extent = renderer.fullExtent()
rect = QgsRectangle(extent)
# print rect.toString()
rect.scale(1.1)
renderer.setExtent(rect)

# set output size
renderer.setOutputSize(img.size(), img.logicalDpiX())

# do the rendering
renderer.render(p)
p.end()

# save image
img.save('test.png', 'png')

QgsApplication.exitQgis()


texts: How to access the different panels in QGIS Composer?

texts:
I accidentally closed the tabs in composer: Composition and Properties, I think it is called.

The question: How can I open them again?


texts: Visualizing business operation time (opening/closing) over course of week in QGIS?

texts:
I'm pretty new to QGIS3.4 and I would like to ask how could I visualize the business operation time (opening/closing) over the course of a week.

I know that I could input an time in the data attributes but
how can I input a data of time over the course of the week with some exception in closing-opening time like the store closes on Sunday,or the there are some different opening-closing time each days.


texts: Reverse gdal_rasterize (-i option) does not work

texts:
When I use:

gdal_rasterize -b 1 -b 2 -b 3 -burn 255 -burn 255 -burn 255 polygon.shp tiff.tiff

I get a new tiff with a white hole in it with the shape of the polygon shp. However, when I use the reverse process:

gdal_rasterize -i -b 1 -b 2 -b 3 -burn 255 -burn 255 -burn 255 polygon.shp tiff.tiff

I get a blank tiff. I expect the old tiff shaped by the polygon (ie. raster values inside the polygon), instead.

Any ideas why this happens? is it possible to clip a raster with a shapefile in this way?


texts: Displaying only the negative or positive contours?(not labels)

texts:
I have a coastal environment and I would like to hide land contours and keep isobaths only(and vice versa probably). One solution would be to superimpose the land polygon, but I do not want to do that. What is the expression I must use in QGIS 3.10?

Additionally dew to very steep bottom slopes isobaths are very dense in some spots and the result is not good. How do I loosen up a little? Any ideas?


texts:
I have a map with a polygon layer containing 100+ features.

I want them to display on the legend and the Layer Tree using the Polygon's ID and then its name i.e: "Id - Name".

I'm using the following steps in QGIS:

Symbology -> Categorized -> concat(Id, ' - ', Name).

However, after following the above steps they are displayed in the Legend and Layer Tree in this order :

1 - Name1
10 - Name10
100 - Name100
101 - Name101
11 - Name11
12 - Name12....
2 - Name2
20 - Name20
3 - Name3

So I read that it's sorting them as strings, and I want to know if there is a way to have them sorted by Id while still display like "ID - Name".


texts: Getting "Failed to open tool" message when trying to run geoprocessing tool from C# code

texts:
I've been experimenting with the Military Aspects of Terrain Template, which has a toolbox associated with it.  I've been using the following WPF Runtime code, which should open a tool in the associated toolbox . . .

private void RunToolButton_Click(object sender, RoutedEventArgs e)
{
        ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Desktop);
        //Set a reference to the IGPCommandHelper interface.
        IGPToolCommandHelper pToolHelper = new GPToolCommandHelperClass() as IGPToolCommandHelper;

        //Set the tool you want to invoke.
        string toolboxName = @"C:\ArcGISForDefense\Intelligence\MilitaryAspectsOfTerrainTemplate\Maps\Toolbox\Military Aspects of Terrain Toolbox.tbx";
        pToolHelper.SetToolByName(toolboxName, "CoverFromFlatTrajectoryWeapons");

        // Generate the array of parameters.
        IArray parameters = new ArrayClass();
        parameters.Add(@"C:\ArcGISForDefense\Intelligence\MilitaryAspectsOfTerrainTemplate\Maps\OperationGDB.gdb\BrigadeAOI");
        parameters.Add(@"C:\ArcGISForDefense\Intelligence\Operational_Environment\Data\MAoT\MAoT.gdb\CombinedVegetationCoverage");
        parameters.Add(@"C:\ArcGISForDefense\Intelligence\MilitaryAspectsOfTerrainTemplate\Maps\OperationGDB.gdb\CoverAreas");

        //Invoke the tool.
        pToolHelper.Invoke(parameters);
}

When I try to run the code, I get the following error message, then it crashes at the "invoke" method . . . "Failed to open tool CoverFromFlatTrajectoryWeapons (Military Aspects of Terrain Toolbox.tbx)"
Am I using the right syntax for opening a tool that is in the "Cover and Concealment" toolset, which is in the "Military Aspects of Terrain Toolbox.tbx" toolbox?  Am I doing anything wrong that would prompt this error message?  I would appreciate any help.

I also wrote the following routine as an alternative, but this also bombs when it gets to the last line of the routine . . .

 public void Run_CoverFromFlatTrajectoryWeapons_Tool()
    {
        // Initialize the geoprocessor.
        //Insert this line before invoking any ArcObjects to bind Engine runtime.
        ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.Engine);
        GeoProcessor gp = new GeoProcessor();

        // Add the toolbox.
        gp.AddToolbox(@"C:\ArcGISForDefense\Intelligence\MilitaryAspectsOfTerrainTemplate\Maps\Toolbox\Military Aspects of Terrain Toolbox.tbx");

        // Generate the array of parameters.
        IVariantArray parameters = new VarArrayClass();
        parameters.Add(@"C:\ArcGISForDefense\Intelligence\MilitaryAspectsOfTerrainTemplate\Maps\OperationGDB.gdb\BrigadeAOI");
        parameters.Add(@"C:\ArcGISForDefense\Intelligence\Operational_Environment\Data\MAoT\MAoT.gdb\CombinedVegetationCoverage");
        parameters.Add(@"C:\ArcGISForDefense\Intelligence\MilitaryAspectsOfTerrainTemplate\Maps\OperationGDB.gdb\CoverAreas");
        // Execute the model tool by name.
        gp.Execute("CoverFromFlatTrajectoryWeapons", parameters, null);
    }

By the way, the above mentioned tool does run well with the specified parameters when running it in ArcMap.

Thanks,
Renee


texts: With the new 2.6 QGIS version I often have problems with raster colors. For some reason certain pixels do not display initially, but if I move the raster a little spatially they show up.


texts: How does iterative clip decide its processing order?

texts:
I have a fishnet grid which I use to iteratively clip several large shapefiles. The output of the clipping tool results in the selected output filename plus an incrementing suffix (_0, _1, _2 ... _n).

I have noticed while testing that the order of clipping is different every time. Once it goes from bottom left to top right, row by row. The next time it is completely random, and the other time it goes from top left to bottom right, column by column.

This has all been on QGIS 2.18, so how does it decide this order? The only factor is the size of the input fishnet grid, but I don't get how this affects the order of clipping.


texts: Field cannot be deleted because it is being used as relationship key?

texts:
I want to delete some fields of featureclasses, but when I run the code  FeatureClass.DeleteField(Field) it works well for some featureclasses but after some time it gives an error such as:

  The field cannot be deleted because it is being used as a relationship key.


texts: Find the CRS of a spatial object from .prj

texts:
I have two spatial objects (coordinates in .shp files) that I want to draw together in R, it is using the same CRS (i.e. to map these using an identical, appropriate, projection)

(The R code to define a projection is similar to:
SpatialPoints(coords, proj4string=CRS('+init=epsg:4326')))

First shape file has an associated .prj that contains:
GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Second shape file has an associated .prj that contains:
PROJCS["WGS_1984_Web_Mercator_Auxiliary_Sphere",GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Mercator_Auxiliary_Sphere"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Standard_Parallel_1",0.0],PARAMETER["Auxiliary_Sphere_Type",0.0],UNIT["Meter",1.0]]

I tried to do a search using key terms in this .prj such as  https://spatialreference.org/ref/?search=WGS_1984 but don't know which of the several projections to pick.

How can I find the proper argument to CRS()?


texts: Label with different fonts and sizes from a field in QGIS 3.8.3

texts:
I received a shapefile from Arcmap.
It has a column for label representation.
Here's an example: 'PP3<FNT name ="SymbolMono BT" size ="11">g</FNT>gc2'.
In Armap the codes are replaced and we have something like the following acronym "PP3ggc2",

with variation in font type and size.

How can I play this label in QGIS 3.8.3?

Is there any way? Else, how may you get around this situation?

Following is an example of the table.


texts: Converting KML to shapefile via ogr2ogr including all attributes

texts:
Here is the example KML file I am using:

<kml xmlns="http://www.opengis.net/kml/2.2">
<Document>
  <name>My Golf Course Example</name>
  <Placemark>
    <name>Club house3</name>
    <visibility>1</visibility>
    <phoneNumber>554545454</phoneNumber>
    <Point>
      <coordinates>-111.956,33.5043</coordinates>
    </Point>
  </Placemark>
</Document>
</kml>

I have named it input3.kml.

I put it, and ogr2ogr in the same folder, then I run this command:

ogr2ogr -f "ESRI Shapefile" output3.shp input3.kml

After that, I open up Quantum GIS and add vector layer, then add output3.shp.

Then I find the output3 layer, right click and open attribute table. How come only two attributes show up? (Name and Description)


texts:
In the GUI there is the Query Builder, but that does not seem to be easily available in Python. The processing (SEXTANTE) module has the qgis:selectbyattribute algorithm, but that only supports selecting by a single attribute. There is also QgsFeatureRequest.setFilterExpression, but apparently that is not available in current stable version (2.0.1 Dufour).

Do I have to write my own?


texts:
If the user moves the mouse pointer to the top left corner of the screen, I want my software to auto-load the three tiles that are closest to the mouse pointer. My assumptions are that the smaller are the map tiles to be loaded, the lesser tiles I will have to load. The question is,

whether there is any existing logic
that I can use for this loading
(something people have used already,
since I need a flexible way to ) and
...
...a flexible way to store the
loaded tiles in C++?. My app is in
C++, and I'm planning to store each
tile as either a class or a struct.
Is there an existing way/a better
way to do it?

EDIT:
Responding to Dan:
I'm still in fact-finding stage, and currently only have a map with roads stored as lines, where each line is composed of a list of line segments specified by a starting lat-long and an ending lat-long. I'm open to using any tiling scheme, but in a simple way, I'm visualizing the map to be divided into a grid of small squares, so the squares/tiles close to the mouse pointer and containing road segments should be loaded into memory. My algorithm would then process the loaded tiles to find the angle of the road and to check if the road forms a junction with other roads etc. My algorithm can't afford to wait for too many loading operations from the hard-disk. The priority is knowing how to store information so that I'll have them in the RAM as quick as possible. Whether to pre-calculate the tile area and store the tiles in a database or whether to extract square portions from the map and cache them in the RAM was a decision I was trying to take when I posted this question.


texts: Transform complex features into simple features GML with HALE

texts:
Background information:

I'm working with complex GML features and would like to publish my data as a WFS service through GeoServer. I know there is already the application schema support (app-schema) extension which provides support for Complex Features in GeoServer WFS.
What I have to do now is to create a mapping file. I'm just trying to figure out what the best options are.

I found out that the HALE Alignment Editor (See here) is able to create a mapping file from simple features through the alignment between the source and the target schemas.

I think the program is very useful to create a mapping file on a visual basis using simple features as input. Nevertheless, I found no information if you can transform your data in the opposite way (Complex into Simple).

Question:

Is there a way to break down the complex features so that I can import them as simple features into the HALE Alignment Editor?

Something like this:
Complex Features -> "Transformation" -> Simple Features -> "HALE Alignment" -> Complex Features + App Schema -> Ready to use with GeoServer


texts: How to avoid very large output rasters using grass algorithms

texts:
I use r.cost.points from the GRASS algorithms in QGIS to create a accumulated travel costs map that is based on a friction surface raster. My problem is that the output file is more than 20x bigger than the input.

The input file is a raster in .tif format and 5MB filesize. The output file is a .tif as well but filesize increases to 110MB.

I don't really understand this increase, because both have the same extend, resolution, Data Type (Float64) and Filetype.

The stored values differ: while the friction map contains only small integer values that inform the travel costs to cross one gridcell (ranging from 1-36), the output map contains the accumulated travel cost values (ranging from 252-46677.400562602). But still both are Float64 so they should be comparable in terms of space requirements AFAIK.

Comparing both Metadata Properties the only difference is the No Data Value which is in case of the input file -1.7e+308  and in case of the output file nan...

My main problem now is that I want to process an input file of 10GB which causes space problems while the output is created... Can anybody explain why this problem occurs and how I can avoid these increases?

You can test this on your computer with this sample data here as the friction surface and here the source layer containing the destinations.

edits
i also tried it with a rasterized version of the source layer but this still creates a File of 98MB.

I use the GUI but the command looks in the console like this.

0   GRASS GIS 7 execution commands
        g.proj -c proj4="+proj=eqc +lat_ts=0 +lat_0=0 +lon_0=-54 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"
        r.external input="/home/user/R/R-projects/3_accessibility_map/output/toy_data/2_friction/friction_1.tif" band=1 output=tmp1480016848937 --overwrite -o
        r.external input="/home/user/R/R-projects/3_accessibility_map/output/toy_data/1_rasterized/towns2.tif" band=1 output=tmp1480016848938 --overwrite -o
        g.region n=-1340309.99994 s=-1451939.99994 e=-555477.497208 w=-667107.497208 res=30.0
        r.cost  input=tmp1480016848937 start_raster=tmp1480016848938 -n max_cost="0" null_cost="0" memory="4000" output=output5ae5c0b98f304e92ada91772b870452b --overwrite
        g.region raster=output5ae5c0b98f304e92ada91772b870452b
        r.out.gdal --overwrite -c createopt="TFW=YES,COMPRESS=LZW" input=output5ae5c0b98f304e92ada91772b870452b output="/home/user/Desktop/test2.tif"

ant the first lines of the output

2016-11-24T20:47:51 0   GRASS GIS 7 execution console output
        Cleaning up temporary files...

        Starting GRASS GIS...

        Executing '/home/user/.qgis2//processing/grass7_batch_job.sh' ...

        WARNING: Datum <unknown> not recognised by GRASS and no parameters found

        Default region was updated to the new projection, but if you have multiple mapsets `g.region -d` should be run in each to update the region from the default


texts: ImportError: DLL load failed: The specified module could not be found. while importing GeoPandas

texts:
How can I fix the error and get GeoPandas to work on my Win10 machine?
Solutions tried:

New environment following instructions here
Uninstalling and reinstalling Geopandas

Detailed error:
import geopandas as gpd
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Program Files\JetBrains\PyCharm 2019.1.3\plugins\python\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
  File "C:\Users\**\Anaconda3\envs\rider_env\lib\site-packages\geopandas\__init__.py", line 1, in <module>
    from geopandas._config import options  # noqa
  File "C:\Program Files\JetBrains\PyCharm 2019.1.3\plugins\python\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
  File "C:\Users\**\Anaconda3\envs\rider_env\lib\site-packages\geopandas\_config.py", line 126, in <module>
    default_value=_default_use_pygeos(),
  File "C:\Users\**\Anaconda3\envs\rider_env\lib\site-packages\geopandas\_config.py", line 112, in _default_use_pygeos
    import geopandas._compat as compat
  File "C:\Program Files\JetBrains\PyCharm 2019.1.3\plugins\python\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
  File "C:\Users\**\Anaconda3\envs\rider_env\lib\site-packages\geopandas\_compat.py", line 8, in <module>
    import pyproj
  File "C:\Program Files\JetBrains\PyCharm 2019.1.3\plugins\python\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
  File "C:\Users\**\Anaconda3\envs\rider_env\lib\site-packages\pyproj\__init__.py", line 49, in <module>
    import pyproj.network
  File "C:\Program Files\JetBrains\PyCharm 2019.1.3\plugins\python\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
  File "C:\Users\**\Anaconda3\envs\rider_env\lib\site-packages\pyproj\network.py", line 10, in <module>
    from pyproj._network import (  # noqa: F401
  File "C:\Program Files\JetBrains\PyCharm 2019.1.3\plugins\python\helpers\pydev\_pydev_bundle\pydev_import_hook.py", line 21, in do_import
    module = self._system_import(name, *args, **kwargs)
ImportError: DLL load failed: The specified module could not be found.


texts: Copyright issues and credits of the base maps provided in the Quickmapservice plugin in QGIS

texts:
I have used QGIS to create maps for a paper that I am going to submit to an academic journal. However, the journal requires that all figures should be under cc-by 4.0 license, while on the QGIS website, it is stated that "If not stated otherwise, all content is licensed under Creative Commons Attribution-ShareAlike 3.0 license (CC BY-SA)".
Does that mean that I cannot use maps created by QGIS for this journal?
Also, I know that every user has to look for copyright permission and how to cite every base map respectively. I have used the ESRI satellite as the base map but found this World Imagery.
Does that mean that I cannot use these maps as the base maps if I am using QGIS instead of ArcGIS?


texts: Connection pooling in OGR PostgreSQL driver

texts:
OGR (GDAL) provides a vector driver for "PostgreSQL" and data in the DB can be accessed with something like the following:
driver = ogr.GetDriverByName("PostgreSQL")
datasource = driver.Open("PG:host=localhost port=5432 user=postgres password=postgres")
layer = datasource.GetLayerByName("geo_table_name")
print(f"{layer.GetName()} has {layer.GetFeatureCount()} feature(s)")

This approach abstracts and obscures the database connection, and as a result we do not appear to have control over connection pooling.
I intend to use OGR PostgreSQL in an API and there could potentially be significant variation in the number of postgres datasource "connections", varying by load. I would like to know if there is an established pattern for managing database connections in OGR, if I am required to implement something myself, or if the idea itself is fundamentally flawed in some way.
Can OGR PostgreSQL connections come from a connection pool, or is there an established approach for ensuring the number of connection attempts does not grow unsustainably under load?


texts: GDAL Application crashes when loading GDAL drivers?

texts:
I'm developing on QT an application that uses GDAL library (which I compiled using MinGW 4.8).
When I use the application on my windows 8 laptop, it works fine, but when I tried using it on several windows 7 computers it crashes as soon as I use the GDAL dll.

Mapper::Mapper() {
GDALAllRegister();
...
}

My program didn't send any error, but when I looked in the eventlog I saw this:

Faulting application name: MainApp.exe, version: 0.0.0.0, time stamp: 0x00000000
Faulting module name: ntdll.dll, version: 6.1.7601.18247, time stamp: 0x521ea8e7
Exception code: 0xc0000374
Fault offset: 0x000ce753
Faulting process id: 0x12fc
Faulting application start time: 0x01cfc5ade29fe9d9
Faulting application path: C:\Users...\release\MainApp.exe
Faulting module path: C:\Windows\SysWOW64\ntdll.dll
Report Id: 3e036714-31a1-11e4-8ac1-bc7737c9ef36
Fault bucket , type 0
Event Name: APPCRASH
Response: Not available
Cab Id: 0
Problem signature:
P1: MainApp.exe
P2: 0.0.0.0
P3: 00000000
P4: StackHash_940e
P5: 6.1.7601.18247
P6: 521ea8e7
P7: c0000374
P8: 000ce753
P9:

Any ideas what may be causing this?


texts: Show content of a file geodatabase on ArcGIS Online

texts: I uploaded a zipped file geodatabase into ArcGIS Online from the content menu. Now I want to see the features inside my geodatabase displayed on the map and I can't find any option to do that.


texts: Why is it that QGIS shows only one attribute table on all layers?

texts: I have checked the "Select Single Feature" before clicking the Attribute Table but even if I choose a different layer, only one attribute table is shown.


texts: GeoServer WMS Layers not overlaying on OSM base-Leaflet

texts:
I would like to understand why I cannot view or see the WMS layers loaded. Based on the investigation I did so far, I am assuming this has to do with base layers projections. My overlays are WGS84, EPSG 4326. Any guidance on how to resolve this?

This is what I currently see:


texts: How to set different label content for different scales?

texts:
Possible Duplicate:
  How to accomplish conditional labeling based on map scale and attribute in QGIS 1.7

I'd like to show different label content at different map scales. Zoomed in, I use the new labeling (ABC-button) to concatenate different values

( "COMPANY_NA"  || '\n' || 'Status Truck: '  ||  "STATUS_TRU"  || '\n' || 'Phone: '  ||   "PHONE_CENT" )

Using a CSV-layer, I could combine with the old label function at the layer properties. Somehow, this does not work as I change the CSV layer to a shapefile, the scale properties of the old label function are not kept.

Unfortunately, as I plan to publish the result on the web client, it is not possible for me to use two layers as in different labels at various scales

Using Version 1.8 on Windows XP


texts: Update for field values from another table field values with condition

texts:
I have two tables in file geodatabase and both are the same structure.

Both tables are having around 140 columns and 4000 rows.

In table 1 there are many fields are "null" and "0" and remaining fields are filled with some values.

In table 2 filled with some values and also few fields with "null" and "0" values

In both table is having "ID" common unique value field available.

Condition is.

I need to update Table 1 based on below scenarios

If fields values available in both then no need update
If table 1 is having "null" and "0" values then update new values if available in table 2.
If table 1 is having some value and if table 2 is having  "null" and "0" then no need to update


texts:
I have the following code to get area of polygon using shapely and pycrs.
I was able to run this code successfully in a jupyter notebook but when I incorporate it into my main code, I had an error show up:

No handlers could be found for logger "shapely.geos"

Seems like I cannot retrieve the area information of the transformed polygon.
Anybody has any idea why?

def getAreaPolygon(polygon, config):
if polygon != 'POLYGON (())':
    wkb = shapely.wkt.loads(polygon)
    #mv = makeValid(wkb)
    try:
        prj_file = config.inputfile.replace(".jpg",".prj")
        print prj_file
    except Exception as e:
        print(str(e))
        print("prj files does not exist")
    fromcrs = pycrs.loader.from_file(prj_file)
    fromcrs_proj4 = fromcrs.to_proj4()
    proj = partial(pyproj.transform, pyproj.Proj(init='epsg:4326'), pyproj.Proj(fromcrs_proj4))

    ##### Code run successful until here.

    projected_area = transform(proj, wkb).area ## error steps
else:
    projected_area = 0.0
return projected_area


texts: Interpolating FROM points TO points (i.e. not to a new raster)?

texts:
Imagine I have two samples of XY points, lets call them A and B, distributed across the same area/extent:

Points A have an attribute, lets call it N. I’d like to estimate (i.e. interpolate) attribute N for points B.

I could simply “copy across” the value from the nearest corresponding points in A to B, but this is no good for the use case I have in mind. Instead I need to interpolate the value, e.g. using IDW or kriging. So far all tools I can think of take point attributes and interpolate them to a new raster (i.e. continuous surface, which I do not want), rather than interpolating to existing points. Does anyone know of any way to accomplish the latter please?


texts: Calculate length of huge vector lines dataset

texts:
I have a MultiLineString vector layer from which I would like to calculate the length.

I have tried to calculate the length with the Field Calculator by putting in the expression $length. The problem is that the dataset is very big (the attribute table consists of about 700,000 features) and my computer or QGIS is not capable of calculating that (at least when using this approach).

Any suggestions of other ways QGIS might be able to calculate the length and deal with the size of the data given the large number of features?


texts:
I'm attempting to import a kml/kmz file generated using Google Tracks into my QGIS 2.4 system (built with GDAL 1.10). KML files currently lose all attributes, and as such I'm attempting to build GDAL with libKML support as I'm led to believe this will solve this issue - as per https://code.google.com/p/libkml/wiki/BuildingAndInstalling
on my Ubuntu 14.04 system.

I've successfully cloned and built the libkml library and tested it as per instructions in the wiki. I'm now trying to build gdal against the libKML library but keep failing.

I'm attempting to build against the latest SVN version as recommended here: http://www.gdal.org/drv_libkml.html

When running configure my system finds the Google KML headers but can't find the libkml library, even when I specify the path where where I built the library. For example on my system libKML was installed at /source/libkml/build

Has anyone had an experience building GDAL with this library on a linux system?

I can't seem to find any specific instructions for Ubuntu though some Fedora instructions have got me this far. As you can probably tell I'm not super confident using linux command line!


texts: Error when using arcpy.ArcSDESQLExecute() to connect to SQL Server - ArcGIS 10.4

texts:
I am able to connect to a MS SQL Server 2012 geodatabase with the arcpy.ArcSDESQLExecute() when using the .sde connection file pre-created in ArcGIS.

I am not able to connect to a MS SQL Server 2012 geodatabase with the arcpy.ArcSDESQLExecute() when not using the .sde connection file pre-created in ArcGIS.

arcpy.ArcSDESQLExecute(server='mymachinename', instance=None, database='mydbname', user='myuser',password='mypass')

It seems as ArcSDESQLExecute can only connect when providing the log-in info to the SDE App Server and not using the Direct Connect:

  RuntimeError: Object: Connection error Entry for SDE instance not
  found in services file

I have also tried to use the standard port for SQL Server (which is 1433) and the default name for SQL Server instance (MSSQLSERVER in my case):

arcpy.ArcSDESQLExecute(server='mymachinename', instance=1433, database='mydbname', user='myuser',password='mypass')

and

arcpy.ArcSDESQLExecute(server='mymachinename', instance='1433', database='mydbname', user='myuser',password='mypass')

This time another error is thrown:

  RuntimeError: Object: Connection error SDE not running on server

I don't have any SDE App Server; I just want to connect to an SQL Server geodatabase with direct connect without first creating .sde connection file.

Is it a limitation of the ArcSDESQLExecute that you cannot connect on-the-fly to a geodatabase or am I missing some parameters here?


texts: Create Intersecting lines and remove dangles

texts:
I have 3 cases shown in the image.

In all of them i want to intersect the lines and remove dangles using python in Qgis python console so afterwards i can make a plugin for this work in qgis.In case 1 i want to extend the vertical line to the horizontal line.In case two i want to extend both of the lines so that they intersect each other.In case three i want to remove the extra portion (remove dangles).

I am able to get the coordinates in the variable.

Can anybody tell me how to do this using python code?


texts: Programmatically styling layer after output is returned from QGIS 3 processing script?

texts:
Hi this may be more of a python question rather than qgis processing script question. I've managed to successfully get the output that I want from my processing script, but I wish to take one step further and style it automatically with a saved qml file.

The code below doesn't work because the return ends my script. How should I go about doing this?

return {self.OUTPUT: dest_id}
iface.activeLayer().loadNamedStyle('Buffer.qml')
iface.activeLayer().triggerRepaint()


texts: Refreshing WMS based on updated CQL Filter in Leaflet

texts:
How do you update a WMS in a LeafletJS Map?

Using a simple text box and submit button, as well as a variable to hold a CQL filter (SQL query text) I've been able to update the CQL filter value on the WMS layer but cannot figure out how to refresh the WMS image.

The text box simply takes a number, and appends it to a text value to create a simple CQL filter (ex: schnum =  450)

The CQL filter works in the Geoserver Openlayers preview, and serveral methods to test have shown it to work as well.

The problem is one of two things happens depending on the methods I've tried:

a new point is added but the old one is not removed
the basemap disappears

I've created a JSFiddle* to hopefully demonstrate this simple example: https://jsfiddle.net/dpsspatial/fqa4y1cd/1/

*NOTE: for some reason the script to update the SQL filter - set_sql() - doesn't seem to be recognized in the JSFiddle though it works locally.


texts: Schema Lock issue when using ITable.AddField via ArcObjects

texts:
During a process I am running, I sometimes need to create a field that is missing in an attribute table. Unfortunately, the program is complaining that the table I am adding the field to has a schema lock from the program itself. So I attempted to add the ISchemaLock method for granting an exclusive lock and I now get the HRESULT error 0x8004022D when attempting to change the lock. Here is what the code looks like.

private void DoWork()
{
        IWorkspace workspace = getWorkspacefromPath(workspacepath);
        IFeatureWorkspace fWorkspace = (IFeatureWorkspace)workspace;
        IFeatureClass fc = fWorkspace.OpenFeatureClass(targetfile);
        if (fc == null)
        {
            MessageBox.Show("No feature class found.", "Error.");
            return;
        }
        ITable table = (ITable)fc;
        int index = table.FindField("PSTATE");
        if (index == -1)
        {
            ISchemaLock schemaLock = (ISchemaLock)table;
            settesttext("No PSTATE field found. Adding field now.");
            try
            {

                schemaLock.ChangeSchemaLock(esriSchemaLock.esriExclusiveSchemaLock);
                table.AddField(fielddata());

            }
            catch (Exception error)
            {
                MessageBox.Show("Field creation error. Closing process." + "\n " + error.Message);
                return;

            }
            finally
            {
                schemaLock.ChangeSchemaLock(esriSchemaLock.esriSharedSchemaLock);
            }
            settesttext("Field creation complete");
            index = table.FindField("PSTATE");
}

How should I go about adding the field to avoid getting a schema lock problem?


texts: Converting Unix epoch time to GMT time using QGIS Field Calculator

texts:
I am trying to convert the 'time' and 'updated' column from epoch to date/time GMT but I am struggling to do so. I know I have to use the field calculator tab but I am a basic at using scripting.
I am using QGIS.


texts: ArcGIS Javascript API - Using Proxy For Specific Layers

texts:
With ESRI's Javascript API - is it possible to make only certain layers go through the proxy?
I know this is possible in the Flex API via proxyUrl property of applicable layers.

Right now i have this code in my initialization function - but all map content is now routed through proxy.

esri.config.defaults.io.proxyUrl = "proxy.ashx";
esri.config.defaults.io.alwaysUseProxy = true;


texts: Writing shapefile with projection defined crashes fiona

texts:
I've just started trying to use Fiona to write out shapefiles. I found the example code below, which seems to work fine

from shapely.geometry import mapping, Polygon
import fiona

# Here's an example Shapely geometry
poly = Polygon([(0, 0), (0, 1), (1, 1), (0, 0)])

# Define a polygon feature geometry with one attribute
schema = {
    'geometry': 'Polygon',
    'properties': {'id': 'int'},
}

# Write a new Shapefile
with fiona.open('my_shp2.shp', 'w', 'ESRI Shapefile', schema) as c:
    ## If there are multiple geometries, put the "for" loop here
    c.write({
        'geometry': mapping(poly),
        'properties': {'id': 123},
    })

However, if I alter this to define a projection for the output shapefile, by changing the open call to:

from fiona.crs import from_epsg
with fiona.open('my_shp2.shp', 'w', 'ESRI Shapefile', schema, crs=from_epsg(3405)) as c:

then Python hangs for a while, and then crashes with the following error:

Python(12985,0x7fff72746180) malloc: *** error for object 0x7fce4aa2c930: pointer being freed was not allocated
*** set a breakpoint in malloc_error_break to debug
[IPythonQtConsoleApp] WARNING | kernel restarted

Does anyone have any idea what's going on here? Am I doing something crazy, is my computer misconfigured somehow, or is this a bug?

I'm running this on OS X 10.8.5, and I've been using GDAL/OGR on this machine for a long time with no problems. The full stack trace, in case that helps, is available here.


texts: Delete subselected features with PyQGIS and processing

texts:
In QGIS 3.10, I have layers in which I select features by location (to find duplicates), using :

selected_features = processing.run('qgis:selectbylocation', {'INPUT':layer1, 'PREDICATE':0, 'INTERSECT':layer2, 'METHOD':0})

'METHOD':0 means I begin a new selection with the selected features

Then, in this selection, I want to select features by expression, using :

subselected_features = processing.run('qgis:selectbyexpression', {'INPUT':layer1, 'EXPRESSION':expression, 'METHOD':3})

'METHOD':3 means I select features in the current selection

Now, I want to delete all the subselected features. I know how to delete features by expression, with the method explained here (Deleting selected features using PyQGIS?) :

with edit(layer1):
    # build a request to filter the features based on an attribute
    request = QgsFeatureRequest().setFilterExpression('"DN" != 3')

    # we don't need attributes or geometry, skip them to minimize overhead.
    # these lines are not strictly required but improve performance
    request.setSubsetOfAttributes([])
    request.setFlags(QgsFeatureRequest.NoGeometry)

    # loop over the features and delete
    for f in layer1.getFeatures(request):
        layer1.deleteFeature(f.id())

But this way selects features in the entire layer1. It does not allow to choose a METHOD as with processing, to make a subselection in the current selection.

Does someone has an idea to delete these subselected features ?


texts:
I've been trying to do project my map in OpenLayers 4 into Massachusetts Mainland (EPSG:26986). I was following this example: http://openlayers.org/en/latest/examples/scaleline-indiana-east.html but after plugging in the MA Mainland values I just get a blank screen. Here is my code:

<script>
proj4.defs('EPSG=26986', 'PROJCS["NAD83 / Massachusetts Mainland",GEOGCS["NAD83",DATUM["North_American_Datum_1983",' +
SPHEROID["GRS 1980",6378137,298.257222101]],PRIMEM["Greenwich",0],' +
'UNIT["degree",0.0174532925199433],PROJECTION["Lambert_Conformal_Conic_2SP"],' +
'PARAMETER["false_easting",200000],' +
'PARAMETER["false_northing",750000],' +
'PARAMETER[""central_meridian",-71.5],' +
'PARAMETER["latitude_of_origin",41],' +
'UNIT["metre"]]');

var map = new ol.Map({
layers: [
new ol.layer.Tile({
source: new ol.source.OSM()
})
],
target: 'map',
view: new ol.View({
projection: 'EPSG=26986',
center: ol.proj.fromLonLat([-71, 42], 'EPSG=26986'),
zoom: 7,
extent: ol.proj.transformExtent([-172.54, 23.81, -47.74, 86.46],
'EPSG:4326', 'EPSG=26986'),
minZoom: 6
})
});

map.addControl(new ol.control.ScaleLine({units: 'us'}));
</script>

Does anyone know why this isn't working?


texts: GDAL ogr2ogr Converting DateTime to String

texts:
I'll preface this by saying I'm brand new to using GDAL. I'm using it in the Linux (Ubuntu) terminal. The GDAL version is 3.0.4.
I'm using ogr2ogr to convert ESRI geodatabase feature classes to GeoJSON files and I've run into a curiosity. Here's the command I ran:
ogr2ogr -f GeoJSON -t_srs EPSG:3310 /mnt/Events/testout1.geojson /mnt/Events/20200910_Event_ArcPro_2_5.gdb.zip Event_Polygon

It gave me this warning: Warning 1: The output driver does not natively support DateTime type for field CreateDate. Misconversion can happen. -mapFieldType can be used to control field type conversion. When I actually inspect the GeoJSON output everything seems fine. The field is of the DateTime type and has the correct values.
But since I'll be doing this iteratively for dozens of files, to be safe I decided to convert any DateTime fields to strings by using the -fieldTypeToString argument. Here's the command:
ogr2ogr -f GeoJSON -t_srs EPSG:3310 -fieldTypeToString DateTime /mnt/Events/testout2.geojson /mnt/Events/20200910_Event_ArcPro_2_5.gdb.zip Event_Polygon

This worked without throwing any warning. But when I inspect the relevant fields in the new GeoJSON file, they're still of the DateTime type! At least that's what QGIS reads it as, and I used ogrinfo to confirm. Anyway, this doesn't actually seem to be a problem since the information is getting carried over properly, but I was curious if anyone had any insight into what's going on here. What if I truly wanted the field to be a string instead of DateTime? If GDAL is leaving the field as DateTime even when I tell it not to, why does it throw a warning about it? I should note that I also tried using -mapFieldType, as the warning suggests, with the same results.