texts: Choosing projection for rasterization of random latitude and longitude data in northern Atlantic?

texts:
I have a dataset in the WGS84 geographic coordinate system that I would like to interpolate along a grid using Python.

It appears gdal_grid is not available through the Python bindings, so I plan to use another 2D interpolation method after rasterizing the data.

To rasterize the data, I plan to re-project it using pyproj, then write it to a raster using the gdal python bindings. The UTM projection coordinate system seems like it would be the best projection, but my data spans the North Atlantic (i.e. multiple UTM zones).

I'm open to suggestion on the interpolation method, but I had planned to use an inverse distance interpolation (as shown here). These are presence absence sightings data from ship cruise tracks, which do not necessarily follow lines of latitude. They span from latitudes in the mid 70's to mid 50's.

What projection would be best to use?


texts: How to create a array variable containing the values of another variable in a table?

texts:
I am new to mapbasic and would like to have  an array variable (sac) in order to access it through the "For variable = 1 To n" function. I would like to create a loop that only takes on the value of a variable (var) in a table and loops through all the values in that Table.

Include "Mapbasic.def"
Declare Sub Main()

Dim sac(100) As Integer
Dim var As Integer

Sub Main()

Sac(1) = 50
Sac (100) = 1200

For sac = 1 to sac(100)

Next

End Sub


texts: Runtime Error: Cannot Open Feature Class in Editor Session

texts:
I am using an update cursor in a script and need to make edits from arcpy.da.Editor.  The error that is being thrown is a Runtime Error where the feature class cannot be opened.

import arcpy
import GetSRConfig
import os

def GetThisDir():
    import inspect
    return os.path.dirname(os.path.abspath(inspect.getfile(inspect.currentframe())))

sDir = GetThisDir()
sSDECon = os.path.join(sDir, GetSRConfig.SDEConnFile)
sdeLayer = os.path.join(sSDECon, GetSRConfig.LN_SDEDATA)
pWorksapce = (sSDECon)

edit = arcpy.da.Editor(pWorksapce)

edit.startEditing(False, False)

srField = "NUMBERCYLA"
vehNameDi = dict ([(sr, vehNam) for sr, vehName
                in arcpy.da.SearchCursor
                (GetSRConfig.LN_SDEDATA, [GetSRConfig.FN_NUMBERCYLA, GetSRConfig.FN_VehiclNam],"""VehicleNam <> ''""")])

with arcpy.da.UpdateCursor(GetSRConfig.LN_SDEDATA, [GetSRConfig.FN_NUMBERCYLA,GetSRConfig.FN_VehiclNam],"""VehicleNam = ''""") as cursor:
    for sr, vehNam in cursor:
        if not sr in vehNameDi:
            continue
        row = (sr, vehNameDi [sr])
        cursor.updateRow (row)

    del cursor
edit.stopEditing(True)


texts: Clip raster by mask layer change the pixels size

texts:
I've a raster that have the pixels of this dimensions:

x 0.2282548198093489467 m
y 0.3594572165075641257 m

I need to clip this raster with a vector because I need a raster with the same extensions of vector.

I use QGIS 3.4 and the process Clip Raster by Mask Layer to do this. At the end of the process I have a raster with this pixels dimensions: 0.345 m.

Why happen this? I expect an output raster with the same pixels dimension of the input raster.


texts: Checking if ModelBuilder model is in use by testing for Schema locks?

texts: Is it possible to get ArcGIS to check if someone else is using a model by checking schema lock and potentially display an error saying that it is already in use and maybe provide information on who it is thats using it?


texts: Can QGIS 3.10.2 w GDAL 3.0 be installed on Ubuntu 18.04 (without being experimental)?

texts:
I just updated QGIS on Ubuntu 18.04 but it didn't update GDAL to 3x, instead it is still at

Running against GDAL/OGR 2.2.3

Is this the correct situation?

Without going into an experimental bleeding edge scenario, can GDAL 3x be put in there?

ta


texts: How much is the error in physical longitude and lagitude (e.g. meters) corresponding to a +-0.0001 error in GPS coordinates?

texts:
If I have a device how give-me a latitude and longitude with 7 decimals (eg. -42.5124478, -22.3006726, longitude and latitude respectively ) and I round it to 4 decimals (eg. -42.5124, -22.3007).

It represents +-0.0001 of possible precision error.

The point is, imagining about a Earth a perfect sphere, I want to know what max limits from rounded decimals in Earth Surface according of rounding (based on this metric answered here Measuring accuracy of latitude and longitude?

how much dislocation of surface area it will represents? and how to calculate it?


texts: Using related layers/tables in background?

texts:
I have postgresql/postgis database setup with one to many relationship. I can successfully set it up in QGIS but i would like to keep related tables hidden from Layers Panel to keep the project "clean" since later on i might be adding more complex databases, value map tables etc..

For example, i have locations layer where each point can relate to several events in events layer. I would like to edit data in form like this so each location can be updated with new events only that events table would not be visible in Layers Panel:

I found solution that works in QGIS 2 but i need to set it up in QGIS 3.

Is there some way to achieve this without going into PyQGIS?

I am just starting with PyQGIS.


texts: QGIS 3.4 Couldn't load SIP module. Python support will be disabled on Windows 10

texts:
I upgraded from QGIS 3.2 to 3.4 and bad things happened.  On start-up I get the error.
Couldn't load SIP module.
Python support will be disabled.

Traceback (most recent call last):
  File "", line 1, in
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python\qgis\__init__.py", line 80, in
    import qgis.gui
  File "C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python\qgis\gui\__init__.py", line 27, in
    from qgis._gui import *
ValueError: PyCapsule_GetPointer called with incorrect name

Python version:
3.7.0 (v3.7.0:1bf9cc5093, Jun 27 2018, 04:59:51) [MSC v.1914 64 bit (AMD64)]

QGIS version:
3.4.2-Madeira 'Madeira', 22034aa070

Python path:
['C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python', 'C:/Users/Cary/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python', 'C:/Users/Cary/AppData/Roaming/QGIS/QGIS3\\profiles\\default/python/plugins', 'C:/PROGRA~1/QGIS3~1.4/apps/qgis/./python/plugins', 'C:\\Program Files\\QGIS 3.4\\bin\\python37.zip', 'C:\\PROGRA~1\\QGIS3~1.4\\apps\\Python37\\DLLs', 'C:\\PROGRA~1\\QGIS3~1.4\\apps\\Python37\\lib', 'C:\\Program Files\\QGIS 3.4\\bin', 'C:\\Users\\Cary\\AppData\\Roaming\\Python\\Python37\\site-packages', 'C:\\PROGRA~1\\QGIS3~1.4\\apps\\Python37', 'C:\\PROGRA~1\\QGIS3~1.4\\apps\\Python37\\lib\\site-packages', 'C:\\PROGRA~1\\QGIS3~1.4\\apps\\Python37\\lib\\site-packages\\win32', 'C:\\PROGRA~1\\QGIS3~1.4\\apps\\Python37\\lib\\site-packages\\win32\\lib', 'C:\\PROGRA~1\\QGIS3~1.4\\apps\\Python37\\lib\\site-packages\\Pythonwin']

I feel that there is an issue with my Python Path but I've tried several things to fix it and I am not there yet.
Does anyone know a fix for this?


texts:
I have some data in a PostgreSQL database, I use pgAdmin4 for the organisation of the data.
I want to create a QGIS project for other users with the data from my database. The users of the project should see the data, load own data to the project, but they should not be able to change my database files.

Therefore I create a user qgis in pgAdmin with this SQL statemant

CREATE ROLE "qgis" with
LOGIN
NOSUPERUSER
NOCREATEDB
NOCREATEROLE
NOINHERIT
NOREPLICATION
CONNECTION LIMIT -1
PASSWORD 'xxxx';

Aftre the creation of the user I give him the privilege SELECT for the different tables in my schema with the statement

GRANT SELECT ON TABLE example_1 TO qgis;

But with this privilege the user can't load the data in QGIS. So I give the user the privilege USAGE. Then the user can load the data in QGIS, but he can also delete data in the database and add new file.
How can I create a user, that can connect to my database, read the data, but nothing more?


texts: Removing distorted edges of TIN surface?

texts:
I'm trying to remove the stretched edges of a TIN surface without having to draw a polygon around the features that I want clipped.

Does anyone know how to do this without free hand drawing?

I was thinking an edge detection method could work.

Here's an example:


texts: OpenLayers3 and GeoJSON data with null features

texts:
I am new to OpenLayers and I'm trying to label features on a map using a GeoJSON data. Here is my code:

var layer = new ol.layer.VectorTile({
    source: new ol.source.VectorTile({
        format: new ol.format.GeoJSON(),
        tileGrid: ol.tilegrid.createXYZ(),
        url: "http://abcd.com/cities/{z}/{x}/{y}.geojson"
    })
})

Sometimes I get "features": null in the GeoJSON data:

{ "type": "FeatureCollection", "features": null }

and features parsing fails with error Uncaught TypeError: Cannot read property 'length' of null. Is it possible to replace null by [] before parsing features?

UPDATE. I found a solution, but got a new issue: I can't add features to the source. Here is my code (I'm using OpenLayers 3.9.0):

var format = new ol.format.GeoJSON();
var source = new ol.source.TileVector({
    tileGrid: ol.tilegrid.createXYZ(),
    tileLoadFunction: function (url) {
        var xhr = new XMLHttpRequest();
        xhr.open("GET", url);
        xhr.onreadystatechange = function () {
            if (XMLHttpRequest.DONE === xhr.readyState && 200 === xhr.status) {
                var json = xhr.responseText.replace("null", "[]");
                var features = format.readFeatures(json, {
                    featureProjection: "EPSG:3857"
                });
                features.forEach(function (f) {
                    console.log(f.getGeometry().getCoordinates()); // features are parsed correctly
                });
                source.addFeatures(features);
            }
        };
        xhr.send();
    },
    url: "http://abcd.com/cities/{z}/{x}/{y}.geojson"
});
source.on("addfeature", function () {
    console.log(source.getFeatures()); // outputs an empty array of the features; it seems that the features haven't been added
});
new ol.Map({
    target: "map",
    layers: [
        new ol.layer.Vector({
            source: source
        })
    ],
    view: new ol.View({
        center: ol.proj.fromLonLat([0, 0]),
        maxZoom: 12,
        minZoom: 1,
        zoom: 6
    })
});

Any suggestions?


texts:
The website http://download.geofabrik.de/north-america.html does not have shapefiles (shp.zip) for sizable states or countries.
How do you convert the available osm.pbf to shp.zip?

I have tried QGIS but it does not export all the shapefiles from (osm.pbf) to (shp).

I have tried R, but again, only exported a fraction of shapes.

library(rgdal)

osmfilename = 'california-latest.osm.pbf'
directory_name = 'shp_directory'
layer_name = 'output_layer'
osm <- readOGR(osmfilename, 'lines')
writeOGR(osm, directory_name, layer_name, driver = 'ESRI Shapefile')

Do you have any other way that I could try?


texts: What is NASADEM and when it will be released?

texts: I heard some time ago about a new global DEM called NASADEM that was going to be released in 2017, and promised to be the best global DEM to date. What exactly is NASADEM and when will it be available?


texts: How do I read and modify heat map value classification on QGIS?

texts: I am trying to create a heat map on human encroachment using the heatmap plugin of QGIS 2.2. I have no problems creating a good looking raster file using a 2 km radius but I can't find my way around the values that are created as a result. I would need values in a scale such as number of points/sq km. Can anybody help me out with this? I know it is a quite basic question but I can't seem to find a solution.


texts: Placing data-defined labels both by expression and manually

texts:
In my project layers (QGIS 3.6.2) I use data-defined label placement. I am trying to set up the placement as follows:

By default, the placement chosen in Label Settings is used (in my case, "around centroid").
If the X and Y fields are populated (i.e., the label has been placed manually), these values are used.
I have a boolean T/F field called "UsePOS" in my data; if it is checked (true), and if the X and Y fields are null, then the X/Y placement uses a calculated point (in my case, a point-on-surface). If it is false, the null X and Y fields are used (i.e., the default placement).

Here's what I tried: in Label Settings, I defined the X placement with the following expression…

if( "X" is not NULL,
    "X",
    if( "UsePOS",
        x(  point_on_surface(  $geometry  )  ),
        "X"
      )
  )

…and the Y field with a like expression. This worked perfectly as expected—by default the label is placed at the centroid; if the box is checked, it jumps to the point-on-surface, and if there is data in the X/Y fields, then the label is pinned there.

However, the problem arises when I try to manually place a label. In that case, because the X/Y placement isn't directly defined by the X/Y fields, QGIS invokes auxiliary storage, which is not what I am using for label placement, and which apparently cannot be suppressed (see my comment in this thread).

So I assume that 1) I need to directly assign label placement to my X/Y fields, 2) I need to use an expression in Field Calculator to populate those fields, and 3) those fields need to be virtual so that they respond dynamically to whether I manually move the label, whether the checkbox is checked, etc. However, entering the above expression in Field Calculator causes QGIS to crash when the label is moved manually, or when I check the checkbox (which I expected, since I don't think a virtual field can reference itself).

The question is, what is the proper way to set up what I want, so that I can still manually move labels, using only data-defined labeling and not auxiliary storage?


texts:
I am displaying icon in feature. I want the size of icon to be dependent to zoom level of map i.e. the icon appears smaller in low zoom levels and bigger in high zoom levels. i used stylemap for that but the icon size does not change on zoom change.

var styleschool = new OpenLayers.Style({
            pointRadius: "${radius}",
            label:"${name}",
            fillOpacity: 0.7,
            fillColor:'blue',
            strokeColor:'blue',
            fontcolor:'blue',
            externalGraphic:'icon/school_color_icon.png'
        },
        {
            context:{
                radius: function(feature){
                    return map.getZoom()*0.9;
                }
            }
        })


texts: ST_3DIntersects says 'no' despite geometries intersecting


texts: Executing script in QGIS after complete load?

texts:
I'm looking for a solution to run a Python script after the complete load of QGIS.
The startup.py is executed during (or before) complete load, so it doesn't work well. I would like to trigger the event "complete load", but I don't know where to look. Where should I look?


texts: How do you aggregate attribute table data at LSOA level to MSOA level?

texts: I have population data stored at the LSOA (Lower layer Super Output Area) but I want to aggregate the data into MSOA (Middle layer Super Output Area). I have the shapefiles for both LSOA and MSOA - each MSOA is broken into 3-5 LSOAs. I have a feeling I need to use "Join attributes by location" but after trying this the population totals for all MSOAs compared to all LSOAs does not add up.


texts: BIOMOD Error message: argument of length 0

texts:
I’m using the package BIOMOD in R to modelling habitat in rodent species. The data formatting works but I'm getting this error message when I run Modelling:

## Conputing the models ##
myBiomodModelOut <- BIOMOD_Modeling(biomodDataL,
                            models = c('GLM','MAXENT','GAM'),
                            models.options = myBiomodOption,
                            NbRunEval = 2,
                            DataSplit = 80,
                            Yweights = NULL,
                            VarImport = 2,
                            models.eval.meth = c('TSS','ROC'),
                            SaveObj = TRUE,
                            rescal.all.models = TRUE)

  ################################################################################ Creating suitable Workdir...

                    ! Weights where defined to rise a 0.5 prevalence !

  -=-=-=-  Lemming  Modeling Summary -=-=-=-=-=-=-=-=-=-

  1  environmental variables ( X2 )

  Number of evaluation repetitions : 3

  Models selected : GLM MAXENT GAM  Total number of model runs : 12

  -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-

  -=-=-=- Run :  Lemming_PA1

  -=-=-=--=-=-=- Lemming_PA1_RUN1    Model=GLM quadratic with no interaction

    Stepwise procedure using AIC criteriaError in 1:ncol(explVar) : argument of length 0

  #

The explanatory variable is an index of productivity (i.e. NDVI). This raster, with explanatory variable (2'x 2' resolution), and vector with the response variable are in the same projection. Do you know why I have this kind of message?

Thanks in advance for your help,
Mary


texts: Combining overlapping polygons to make a raster

texts:
I'm working in Arcmap10.3. I have some polygons layers that over lap that I want to merge to one polygon layer.  Then I want to convert it to a raster in the end.  I want the top layer of the 4 input layers to not be disturbed or converted to the same class as the lower layers.  How do I go about doing this?

I've tried using the merge tool for the 4 polygon shapefiles. But it just runs into infinity.


texts: I am doing a project that is looking at NDVI in certain areas of Africa in Google Earth Engine. I want to take the areas where I have made an NDVI band and then calculate Moran's I so that I can export it to my datasheet. I know how to export it, but I have found no documentation on how to calculate Moran's I within Google Earth Engine. Can it be done? Or will I have to export the data out to R or Python and calculate it there? If so, how do I go about doing that?


texts: Exporting to DXF with properties QGIS 3.12

texts: I am using QGIS 3.12. Is there a way I can export my project to DXF including the information in the properties table, so I will be able to see the properties of the feature in AutoCAD?


texts: Performing gap analysis with line and polygon layer in QGIS

texts:
How do I perform gap analysis with line and polygon layer in QGIS?
I have 2 vector layers. The first vector layer has line features and second layer has polygon features as in the following images.

I am trying to perform a gap analysis that is cutting the lines at segments the polygons do not overlay the lines in QGIS.
Shapefiles of line and polygon layers are accessible from the following link: https://drive.google.com/drive/folders/11GhQODoPby9SCnp9WcydOBJHVS3O1Lsv?usp=sharing


texts: Using primary key in ArcSDE Enterprise Geodatabase with MS SQL Server?

texts:
Does setting a primary key constraint in a field obstruct the compression of a versioned geodatabase?

I read someone mentioning that the primary key in a versioned geo-database will cause problem while compressing the geo-database.

I have a cadastral database in ArcSDE 9.x (MS SQL Server 2008). I have a field parcel id which  is the primary key.


texts: Get the intersection of a polygon with a generated buffer in postgis

texts:
I have two tables. One has points in EPSG:4326 and the other one has polygons also in EPSG:4326.
I loaded these two as layers in geoserver and show them using openlayers (as shown in image)

Now I want to write a query in postgis in which I find the intersection area of each green point (buffered by 200 meters) with the polygons.

I have written this query:

SELECT c.name, d.name FROM points c, polygons d WHERE ST_INTERSECTS(d.geom, ST_BUFFER(c.geom,500))

But I get no results. If I increase the buffer size to 5000000 then I get results.

So the question is if there is a better way to do this without using a Buffer.
If there is projection issue in which I am missing.

EDIT

I followed the advice of @HimBromBeere and transformed my geometries to EPSG: 2100 (which is in metres).
The problem still remains as executing this query I get as a result 0, which doesnt make any sense.

SELECT ST_AREA(
      ST_INTERSECTION(ST_Buffer
      (ST_Transform(ST_SetSRID(c.geom, 4326),2100),200),
      ST_Transform(ST_SetSRID(d.geom, 4326),2100)))
      FROM categoriesdata c, dimoi d


texts: Making nonoverlapping service area polygons using QNEAT3; service area polygon errors

texts:
I am doing a project in which I am determining which areas of New York City can be accessed within a five- or ten-minute walk of a subway entrance. I have done this in the past in ArcGIS, and am trying it out in QGIS. The goal of the project is to determine which areas of the city would become accessible within a five- or ten-minute walk if closed subway entrances, of which I have accumulated a list here, were reopened. I used the QNEAT3 tool to first look at the service areas for open subway entrances, based on the city's street network.
In ArcGIS, I could create non-overlapping service area polygons to allocate areas to the closest subway entrance. I have not seen any way of doing this in QGIS. When using QNEAT3, all the areas within, for example, five minutes of an entrance, are in a single layer, which is not very helpful, as can be seen in this screenshot below.

For some reason, the polygons for areas within a closer range of an entrance, i.e. 60-120 seconds, are sometimes in the incorrect spot, being underneath layers like the 300-480 second polygon layer. If you compare the two images below, you will see that if I uncheck this layer, the 60-120 second layer is visible. In ArcGIS, I had no issues with this at all.
How do I remedy this?

Another issue is that there are spikes in the polygons, especially those in the 480-540 and 540-600 second layers, which should not be there. These areas are clearly not within ten minutes of an open subway entrance. I have no idea what error is causing them.
Does anyone have an idea?
Here is an overview showing multiple errors:

This major one occurs over a body of water:

This one for some reason has a polygon connecting Lower Manhattan with Brooklyn over New York Harbor.

You also have a big one here in Queens.


texts: Openlayers 2 - how to effectively debug rendering failed

texts:
When rendering multiple layers on a map, only one layer renders, the last layer in the array when looking at myMap.layers

On the layers which don't render, the features which are not rendered are held in myMap.layers.unrenderedFeatures

How can I debug WHY they weren't rendered? Is there any more information anywhere that will highlight the issue?

I haven't included code at this stage as it's a messy legacy system and would like to debug myself if possible to get a better understanding of what it's doing and why it's not rendering.

I did look at the docs but they're not particularly helpful


texts: Automating importing of symbology into feature classes in ArcMap?

texts:
I'm using ArcMap 10.

I'm importing several shp files with the same schema using the FeatureClassToGeodatabase tool in Arcpy and during the process I would also like to apply some standard symbology to the feature class.  I was thinking I could use the ApplySymbologyFromLayer but it always errors out when I try to do it with a feature class in a GDB.  It works fine with shp files when using the ArcMap built-in tool.  So something is different when using a feature class within a GDB for the ApplySymbologyFromLayer.

So can I not keep a symbology stored with feature class when it's sitting in the Geodatabase?  Meaning the symbology has to applied when I add the feature class to the data frame?  Or can I use a Representation by using the Add Representation tool?  When I use this it seems to always import the same symbology classification values from the original feature class the representation was created with.


texts: Change in PostGIS database password prompts for user credentials on QGIS start

texts:
I am using QGIS 3.4.6 and I have a project file which was set up with a PostGIS database connection with basic authentication prompting for a username and password. This was stored as plain text in the project file.

Recently the person who manages the PostGIS database changed the password but kept the username the same. This has prompted QGIS to ask for the database credentials whenever the project file is now opened or any past projects that I have worked on.

If I enter the credentials that a connection is made so I know that this is all correct, however if I edit the database connection with the new password and re-save the project file I am still asked for the credentials on start-up.

I have saved the new credentials and stored them as plain text in the project file as before. Interestingly, when I open the project and check the credentials the new ones have been saved.

Does anyone know how I can overcome this and get my project file to remember my user credentials?


texts: How to create a drop-down list of geoprocessing tasks?

texts:
I need to construct a drop down menu with each choice enabling a different geoprocessing task to run over my script.

<select id="myList" onchange="environment()">
   <option>Select..</option>
   <option value="http://********/arcgis/rest/services/EDDM/selectZIP/GPServer">DEV</option>
   <option value="http://********/arcgis/rest/services/EDDM/selectZIP/GPServer">SIT</option>
   <option value="http://********/arcgis/rest/services/EDDM/selectZIP/GPServer">CAT</option>
   <option value="http://********/arcgis/rest/services/EDDM/selectZIP/GPServer">PROD</option>
</select>

They each run the same process albeit at different levels of production. I have this code that is a placeholder and uses just one value.

function selectZIP(Zip_code){
    gp = new esri.tasks.Geoprocessor("http://********/arcgis/rest/services/EDDM/selectZIP/GPServer/selectZIP");

How can I make a drop-down that will allow the user to select which GPServer they need and use it on the same script?


texts: Simple solution for escaping the QGIS locked into fullscreen mode on Ubuntu 12.04

texts:
"

No crash in this case, but I made the mistake of toggling to fullscreen mode on my netbook running Ubuntu 12.04. It's now impossible to escape this "No Menus" view, and it seems to have disabled the popout unity launcher as well. Steps I've already taken:

Alt-Tab to move to a different program, then back - QGIS remains in fullscreen
Alt-V, Ctrl-F and F11 get no response
I've unistalled QGIS, deleted the /.qgis/ profile data folder, and reinstalled. Still launches fullscreen

To add to the weirdness, it now launches - still in fullscreen - into a different workspace than the active one.

I realize this has the sound of a Unity problem, but it's just this program that has a malfunctioning UI.

Thanks to all the crew; this in no way changes my opinion of QGIS as an essential tool.
qgis ubuntu
shareedit

asked Sep 13 '12 at 17:28
Bill Morris
1179
"

***I got the same problem mentioned above and also found a simple solution. I could't post any answer due to lack of reputation in this forum (I am a newbie here).
The solution for escaping the Full Screen trap of QGIS on Ubuntu 12.04 is

(1) - press Alt + F4

(2) - QGIS will ask you about Close without saving, cancel or save, however, when you move your mouse arrow to the top of the screen the menu will appear again, then you go to "view" and click the last option "toggle full screen".

That's all!*


texts:
I have found some info on here regarding opening a KML file (with groundoverlay) from a url in Openlayers. In order to open a KML with groundoverlay, it must first be parsed. However, after I GET, I am getting back a downloadable KML file. From the response in FireBug, I just see a bunch of giberish as the response! I want the KML(XML) so I can get the info out of it. Does anyone know exactly how to parse a downloadable kml from a url? Thanks!

Here is my code below. The url takes me to a kml that I can see in Google Earth. However, the response I get back is just a bunch of symbols according to Firebug, when I send as a GET. Any idea what I am doing wrong here to get the response in xml format?

      kmlURL='http://capstonedev/cgi-bin/proxy.cgi?url=http%3A//10.60.26.64%3A8080/viewshed.kml%3FREQUEST%3DExecute%26DataSource%3DSRTM3%26Operation%3DLegionViewshedOperation%26Parameters%3DobserverCoord%3A-113+42%2CouterRadius%3A40000%26FORMAT%3DKML';
OpenLayers.Request.GET({
    url: kmlURL
    ,
    callback: function(request) {
        if(request.status==200)
           alert(kmlGroundOverlayParse(request.responseText));
        else
        {   alert(request.status + "::" + request.responseText);
            console.log(request.status + "::" + request.responseText);
        }
    }
})

    function kmlGroundOverlayParse(strKML)
    {
    var NL = "\r\n";
    var parserXML = new OpenLayers.Format.XML();
    try {
    objKMLxml = parserXML.read(strKML);
    } catch(err) {
    console.log(err.message);
    //alert(err.message);
    return;
    }

    var Msg="";
    var objLayerList = parserXML.getElementsByTagNameNS(objKMLxml, "*", "GroundOverlay");

    for(var i=0;i<objLayerList.length;i++)
    {

    var derName = parserXML.getElementsByTagNameNS(objLayerList[i], "*", "name");

    var dasIcon = parserXML.getElementsByTagNameNS(objLayerList[i], "*", "Icon");
    var dieHREF = parserXML.getElementsByTagNameNS(dasIcon[0], "*", "href");

    var dieBox = parserXML.getElementsByTagNameNS(objLayerList[i], "*", "LatLonBox");
    var derNorden = parserXML.getElementsByTagNameNS(dieBox [0], "*", "north");
    var derWesten = parserXML.getElementsByTagNameNS(dieBox [0], "*", "west");
    var derSueden = parserXML.getElementsByTagNameNS(dieBox [0], "*", "south");
    var derOsten  = parserXML.getElementsByTagNameNS(dieBox [0], "*", "east");

    Msg += "name:" + derName[0].firstChild.nodeValue + ", ";
    Msg += "url:" + dieHREF[0].firstChild.nodeValue.substr(0,20) + ", ";
    Msg += "BOX:" + derWesten[0].firstChild.nodeValue + ", ";
    Msg += derSueden[0].firstChild.nodeValue + ", ";
    Msg += derOsten[0].firstChild.nodeValue + ", ";
    Msg += derNorden[0].firstChild.nodeValue + ", ";

    Msg += NL;

    var title = derName[0].firstChild.nodeValue;
    var url = dieHREF[0].firstChild.nodeValue;
    var projection = "EPSG:4326";

    var ptmin = [derWesten[0].firstChild.nodeValue,derSueden[0].firstChild.nodeValue];
    var ptmax = [derOsten[0].firstChild.nodeValue,derNorden[0].firstChild.nodeValue];

    var bounds = ptmin[0] + "," + ptmin[1] + "," + ptmax[0] + "," + ptmax[1];
    var size = "100,100"; //uninteressant, notwendig für addLayer??

    addKMLGroundOverlay(i, title, url, projection, bounds, size);
}

Msg += objLayerList.length;
parseMsg = Msg;

}


texts:
I have a code snippet that I can't figure out:

'Define a calculator for NewArea field
    pCalculator = New Calculator
    With pCalculator
        .Cursor = pCursor
        .PreExpression = "Dim dbarea as double" & vbCrLf _
                       & "Dim pArea as IArea" & vbCrLf _
                       & "Set pArea = [Shape]" & vbCrLf _
                       & "dbArea = pArea.area"
        .Expression = "dbArea"
        .Field = "NewArea"
    End With

    'Calculate the field values of NewArea field
    pCalculator.Calculate()

What is pCalculator calculating here?  What is [Shape] in this context?  I'm migrating this old VBA code to Python, and I'm stuck at this segment.

A more general question, what is the significance of the square brackets in the pre-expression?  Are they placeholders for the expression?


texts: VTK format display problem in Paraview+GRASS

texts: I would like to create geological 3D model from lithological layers of an area.I have converted the raster file to VTK(r3.to.vtk command) in GRASS and intended to display it in Paraview. Although I open the vtk file and click apply button, the program(Paraview) always freezes. Why does it happen?


texts: Why did not I manage to reproject layer?

texts:
I have a global land cover map, which has the following properties: Name    ESACCI-LC-L4-LCCS-Map-300m-P5Y-2010-v1.6.1
Path    C:\QGIS\ESACCI-LC-L4-LCCS-Map-300m-P5Y-2010-v1.6.1.tif\ESACCI-LC-L4-LCCS-Map-300m-P5Y-2010-v1.6.1.tif

CRS EPSG:4326 - WGS 84 - Geographic
Extent  -180.0000003392726740,-89.9999950369175963 : 179.9999897345625186,90.0000000000000000
Unit    degrees
Width   129600
Height  64800
Data type   Byte - Eight bit unsigned integer
GDAL Driver Description GTiff
GDAL Driver Metadata    GeoTIFF
Dataset Description C:\QGIS\ESACCI-LC-L4-LCCS-Map-300m-P5Y-2010-v1.6.1.tif\ESACCI-LC-L4-LCCS-Map-300m-P5Y-2010-v1.6.1.tif
Compression LZW
Band 1  •   STATISTICS_APPROXIMATE=YES
•   STATISTICS_MAXIMUM=220
•   STATISTICS_MEAN=188.20256246321
•   STATISTICS_MINIMUM=10
•   STATISTICS_STDDEV=53.693285876722
•   STATISTICS_VALID_PERCENT=100
More information    •   AREA_OR_POINT=Area
•   Authors=Pierre Defourny ,Eric Van Bogaert
•   Copyright=Copyright ©UCL Geomatics, BELGIUM 1999-2012
•   Generate by=gdal_ocl
•   process begin time=2015-08-24T20:26:40
•   process finish time=2015-08-24T20:51:20
Dimensions  X: 129600 Y: 64800 Bands: 1
Origin  -180,90
Pixel Size  0.002777777701186999893,-0.002777777701186999893

To clip continents I used shp file of countries, which has following properties: Fixed geometries
Source

MultiPolygon?crs=EPSG:4326&field=OBJECTID:double(19,0)&field=COUNTRY_NA:string(254,0)&field=Shape_Leng:double(32,15)&field=Shape_Le_1:double(32,15)&field=Shape_Area:double(32,15)&uid={afd52bdd-db05-480d-acf5-932252071748}

Storage
    Memory storage
    Comment

Encoding
UTF-8
Geometry
Polygon (MultiPolygon)
CRS
EPSG:4326 - WGS 84 - Geographic
Extent
-179.9999999999999716,-89.9999999999999432 : 180.0000000000000000,83.6333923340000638
Unit: degrees
Feature count: 279

I collected African countries to be used as mask. Africa was clipped and the new layer has the following properties:

Path   C:\Users\User\AppData\Local\Temp\processing_049cf2b8a5b041639ace295d17d3d405\1c8919fd7e01488e8d7dacc10cd06c93\OUTPUT.tif

CRS
EPSG:4326 - WGS 84 - Geographic
Extent
-25.3611104118373092,-46.9805542601757296 : 73.7444424111124590,37.5444434092434918
Unit
degrees
Width
35678
Height
30429
Data type
Byte - Eight bit unsigned integer
GDAL Driver Description
GTiff
GDAL Driver Metadata
GeoTIFF
Dataset Description
C:/Users/User/AppData/Local/Temp/processing_049cf2b8a5b041639ace295d17d3d405/1c8919fd7e01488e8d7dacc10cd06c93/OUTPUT.tif
Compression

Band 1
STATISTICS_APPROXIMATE=YES
STATISTICS_MAXIMUM=220
STATISTICS_MEAN=115.01767969776
STATISTICS_MINIMUM=10
STATISTICS_STDDEV=74.622799277297
STATISTICS_VALID_PERCENT=30.52
Band 2

More information
AREA_OR_POINT=Area
Authors=Pierre Defourny ,Eric Van Bogaert
Copyright=Copyright ©UCL Geomatics, BELGIUM 1999-2012
Generate by=gdal_ocl
process begin time=2015-08-24T20:26:40
process finish time=2015-08-24T20:51:20
Dimensions
X: 35678 Y: 30429 Bands: 2
Origin
-25.3611,37.5444
Pixel Size
0.002777777701186999893,-0.002777777701186999893

To measure area I used r.report plugin. The result is the following:
I reprojected the clipped layer to 102022 and r.report plugin was running again. The result is the following:
What mistakes have I make?


texts: Google Earth Engine: Accessing image properties inside .map() function

texts:
Background
I would like to access properties of images within an ImageCollection, using the .map() function. Here is an example for an image collection, using yearly rainfall for Tanzania (link to code):
// Define ROI
var TZ_outline = ee.FeatureCollection('USDOS/LSIB_SIMPLE/2017')
    .filterMetadata('country_na', 'equals', 'Tanzania');

// Define time range
var startyear = 2015;
var endyear = 2019;

// List years
var years = ee.List.sequence(startyear,endyear);

//Aggregate total rainfall over year
var yearTotal = ee.ImageCollection.fromImages(
    years.map(function(y) {
        var start = ee.String(ee.Number(y).subtract(1).int()).cat('-05-01') // start rain year from May April-1
        var end = ee.String(ee.Number(y).int()).cat('-04-30')

        return ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD')
                  .map(function(im){return(im.clip(TZ_outline))})
                  .filterDate(start, end)
                  .sum()
                  .set('year', y);
    })
    .flatten()
);

For a single image, I would use .get() to access a property, and this works for my example:
print(yearTotal.first().get('year'))

Using this code within the .map() function returns an error:
var col = yearTotal.map(function(img) {
 var year = img.get('year')
 print(year)
 return(year)
});

Question
How do I make this work? Does it have something to do with server vs. client side?


texts: Calculating one cell in an attribute table in ModelBuilder

texts:
There are two point feature classes--let's call them Established_points and Newcomer. Established_points has, among others, two attribute fields called (for example) ALL_SIGHTINGS and CONFIRMED_SIGHTINGS. In ModelBuilder, I need to append Newcomer into Established_points. Because Newcomer doesn't have the two fields mentioned above, the values in those fields for that record in Established_points will be Null after the appending. I need to calculate the values for those fields, however (based on something else in my model).

What is the best way to go about this? After appending the two feature classes, should I use some kind of cursor, perhaps, to go through the ALL_SIGHTINGS and CONFIRMED_SIGHTINGS fields and calculate any cells that are Null? (There should only be that one cell for each field.) Or is there some easier way, or some workaround with which to avoid the issue altogether?


texts: How to handle the New ESRI Ocean Basemap

texts:
Some of you may have seen this recent (4/28/14) ESRI blogpost - The New Ocean Basemap.

ESRI, in their infinite wisdom, has split up their Ocean Basemap into 2 separate services.  One with just the map tiles.  And one with just the labels.  And they will 'retire' the old basemap in 2015.

My first thought:  Cool, now I'll be able to drop the labels from their basemap if I want.

My second thought:  But what about those times when I don't want to drop the labels...?

With basemaps being mutually exclusive, am I really going to need to build special circumstances into my mapping applications (to say whenever the ESRI Ocean Basemap is active, make sure this other ESRI Ocean Label overlay/layer is also active...etc).

What a PITA.

Are they really making it this difficult to display their Ocean Basemap with the Labels???

Have I mis-read their blogpost?
Am I missing the better/obvious solution?

BTW I've created an openlayers example in a jsfiddle so you can see the 3 ESRI offerings in place (Original Ocean Basemap, New Ocean Basemap, and the New Ocean Labels).


texts:
I am trying to find the number of pubs for each London LSOA.

I do this with London's LSOA boundaries shapefile which I imported as a SpatialPolygonsDataFrame and a SpatialPointsDataFrame of pubs taken from this website with the following query:

/*
This has been generated by the overpass-turbo wizard.
The original search was:
“pub in london”
*/
[out:json][timeout:30];
// fetch area “london” to search in
{{geocodeArea:london}}->.searchArea;
// gather results
(
  // query part for: “pub”
  node["amenity"="pub"](area.searchArea);
  way["amenity"="pub"](area.searchArea);
  relation["amenity"="pub"](area.searchArea);
);
// print results
out body;
>;
out skel qt;

when I run:

result  <- over(pubs, LondonLSOA)

Nothing happens. I think I am supposed to get a new column in results but I dont see that.

Can someone advise?


texts:
How to download NOAA AVHRR 1 km resolution sea surface temperature data?

I need a SST 1 km image of the coast of Netherlands with at least 20 % valid pixels.

After searching and searching I find the NEODAAS portal: http://www.neodaas.ac.uk

I find this image: http://www.neodaas.ac.uk/data/browse/file_info.php?image=AVHRR/uk/2017/04/09/09apr170324uksstp.gif

However I cannot figure out how to actually download the image in SST values? Any suggestions? Or maybe a suggestion for another way of acquiring the image?


texts:
HUC12_intersects is a SpatialPolygonsDataFrame object with over 5000 polygons

library(rgdal)
library(sp)
library(rgeos)
library(raster)
library(maptools)

HUC12_intersects@data
AREA_ACR_1 GEOL_CODE new_soil HRU_ID
22093.3         2        2     10
22093.3         2        2     12
22093.3         2        2     10
22093.3         2        2     10
 1033.8         1        1      3
22093.3         2        2     10
22093.3         2        2      9
22093.3         2        2     10
22093.3         2        2     10
22093.3         2        2     10

HUC12_intersects@polygons[[1]]

An object of class "Polygons"
Slot "Polygons":
[[1]]
An object of class "Polygon"
Slot "labpt":
[1] 2007122.2  840584.2

Slot "area":
[1] 12538.28

Slot "hole":
[1] FALSE

Slot "ringDir":
[1] 1

Slot "coords":
         x        y
 [1,] 2007095 840527.3
 [2,] 2007049 840514.4
 [3,] 2007028 840547.0
 [4,] 2007061 840587.5
 [5,] 2007132 840631.9
 [6,] 2007174 840658.8
 [7,] 2007195 840639.8
 [8,] 2007207 840625.4
 [9,] 2007208 840594.1
 [10,] 2007151 840559.9
 [11,] 2007095 840527.3

 Slot "plotOrder":
 [1] 1

 Slot "labpt":
 [1] 2007122.2  840584.2

 Slot "ID":
 [1] "109"

 Slot "area":
 [1] 12538.28

The output from HUC12_dissolve is out of order

HUC12_dissolve<-aggregate(HUC12_intersects,by="HRU_ID")

HUC12_dissolve@data
  HRU_ID
1       1
2      10
3     100
4      12
5      14
6      15
7       2
8     202
9       3
10      5
11      7
12      8
13      9

I've tried sorting the data.frame and sorting the polygons to match the sorted data.frame but this does not preserve the spatial aspect of the polygons

order_HRUs<-order(HUC12_dissolve@data[,1],decreasing=F)
HUC12_dissolve2<-HUC12_dissolve
cnt<-0
for(i in seq(1,length(order_HRUs))){
  cnt<-cnt+1
  HUC12_dissolve2@polygons[[cnt]]<-HUC12_dissolve2@polygons[[which(order_HRUs==cnt)]]
}

This creates polygons with gaps and they are not labeled with the correct "HRU_ID"

I want to sort the data.frame, and I want the sorted data.frame "HRU_ID" to match the polygons


texts: How to change the coordinate scale of SpatialPolygons in R?

texts:
I hope to plot SpatialPolygons on satellite raster images. The coordinates of SpatialPolygons are converted to be in the range of x = [0,1] and y = [-1,0]. I know that this overlay plot can be done in python. I am curious if we can do this in R as well?

Upate:

> print(readWKT(wkt))
class       : SpatialPolygons
features    : 1
extent      : 0, 0.009158, -0.009045, 0  (xmin, xmax, ymin, ymax)
coord. ref. : NA

> print(raster("three_band/6040_2_2.tif"))
class       : RasterLayer
band        : 1  (of  3  bands)
dimensions  : 3349, 3391, 11356459  (nrow, ncol, ncell)
resolution  : 1, 1  (x, y)
extent      : 0, 3391, 0, 3349  (xmin, xmax, ymin, ymax)
coord. ref. : NA
data source : three_band\6040_2_2.tif
names       : X6040_2_2
values      : 0, 65535  (min, max)


texts:
I have used the distance matrix tool for a large point set (100 000+ points) to find three closest POI's from another layer and their distances. The distance matrix outputs a table much like this:
source   target_id   distance
1        2          100
1        3          110
1        4          120
2        1          100
2        3          200
2        4          210

My point layer is now even bigger because each original point has three separate features which only difference is the target_id and the distance from source to target. I would like to have the target_id and the distance in separate fields in one feature per source_id.
source   target_1   distance_1   target_2   distance_2   target_3   distance_3
1        2          100          3          110          4          120
2        1          100          3          200          4          210
...

I have tried all the basic tools I can use in QGIS without any luck. Seems like I am trying to do a "many-to-one" join, which can not be made with join-by-attribute. Merging or dissolving the layer to itself doesn't help either. I suspect this could be a simple SQL join operation but I have no idea how I could perform it.


texts: How to symbolize points using line symbols of varying rotation and width?

texts:
I've been along a stretch of coast and recorded the orientation and thickness and locations of several igneous intrusions. I've managed to set up a project in QGIS with a point for each intrusion, each point is represented by a line which is rotated so I have the position and orientation.

What I'm trying to do change the thickness of each line to the thickness I recorded and form a sort of pseudo-polygon rectangle (the length of lines is just an arbitrary value) but I think my problem is I'm using rotated points represented by a line symbol :/


texts: How to create a symbol with directional arrows in QGIS 3?

texts:
This is the same as this question for ArcGIS: Creating Directional Flow Arrows for sewer lines? but for QGIS.

I tried to follow the same approach as in ArcGIS, but I keep missing the mark.  I think I am close with a Marker Line.


texts: Not able to export Georeferenced PDF from ArcGIS 10.5.1

texts:
I am trying to export a geo-referenced pdf from ArcGIS 10.5.1.

I have the checkbox ticked in the advanced export settings (export map geo-reference information), but when I open the map using Adobe Acrobat Reader and try to use the measuring tools no georeferenced information is available i.e. The Lat and Long fields are blank.

The same issue has occurred across a variety of maps I have made with varying coordinate systems that are properly functioning within the MXD document itself.


texts: Create vector layers in sql in spatialite db with PyQGIS / QtSql

texts:
I have a vector layer in a spatialite DB structured as follows:

In PyQGIS I try to create a derived vector layer from this like

from PyQt4.QtSql import *
db = QSqlDatabase.addDatabase('QSPATIALITE')
db.setDatabaseName('C:\\tmp\\test.sqlite')
db.open()
q = QSqlQuery('create table nat2000_hetzles_11 as select * from nat2000_hetzles where maes_2_06 = 11', db)
q._exec()

(note that the query is just a simple example, as a matter of fact it might consist of far more complicated conditions, joins, etc.)

The table is created but I'm facing some issues. GEOMETRY column is empty and QGIS just recognizes the result as a table (not vector layer).

The latter can be fixed by registering the geometry column in geometry_columns table like

insert into geometry_columns values ('nat2000_hetzles_11', 'geometry', 3, 2, 31468, 0)

The GEOMETRY seems to need special treatmant, but I could not figure out what.

How is the GEOMETRY created in statements like create my_vec as select * from my_source_vec where...?


texts: How to horizontally align irregularly shaped point symbol to vector layer in QGIS?

texts: Fairly new to QGIS and searching a method to graphically align my point symbols in the direction of the vector. The symbol represents a water valve in an irregular shape in close proximity to a vector water network--is there a way to align each symbol adjacent/ horizontal to a water network vector layer?


texts: Is "cell" in IDRISI (for ex. in area calculation) the same thing as pixel? Or what else it is?


texts: How to calculate area of shapefile with a merged layer

texts:
so i have 'save as''d the new project and layer with my new Custom CRS project string being '+proj=laea +lat_0=70 +lon_0=0 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs'.

However when I calculate area using the field calculator the numbers still don't make sense. Some of them are even in the negative. I have included a screen shot of the layer I am trying to calculate area for. Please note the column titled 'Area_Km2' are the original areas and the column titled 'Area m2' are the new areas which I am trying to calculate.

FYI I am working with a merged layer and I only clipped area in one row namely that titled 'Va2', It should be the only one with changed values. I am assuming the area will be in meters from the definition string.


texts: Using GEOMTRANSFORM twice on class in MapServer?

texts:
I want to know if and how it is possible to use GEOTRANSFORM two times on a class.
1. Creating a bounding box "bbox" of an object
2. Getting the "centroid" of that bounding box to display a vector symbol

Some of the documentation for MapServer GEOTRANSFORM is listed here, but i haven't been able to put the two functions together.

Maybe I should use this example GEOMTRANSFORM (simplify(buffer([shape], 20),10)) and translate it into something like GEOMTRANSFORM (centroid(GEOMTRANSFORM([bbox]))) . I have tried many combinations but so far none have worked. The layer settings in which I want to implement the transformation looks something like this:

LAYER
 GROUP "23Vandforsyning"
 NAME "Ventiler"
 TYPE LINE
 STATUS OFF
 CONNECTIONTYPE OGR
 CONNECTION "DanVand\DVG_Ventil.tab"
 MINSCALEDENOM 0
 MAXSCALEDENOM 3000
 FILTER (("[Status]" != "Ikke i brug") AND ("[Status]" != "Sløjfet") AND ("   [Status]" != "Planlagt"))

CLASS
    EXPRESSION ("[KompFunktion]" == "Stophane")
    STYLE
        GEOMTRANSFORM "centroid"
        ANGLE [VinkelSymbol]
        OUTLINECOLOR 0 0 0
        SYMBOL 'Firkant_m_kryds'
        SIZE 6
        MINSIZE 6
        MAXSIZE 6
    END
END

END


texts: How to set the default values for the ScaleBar in QGIS?

texts:
I often insert scalebars (that's quite common) in my maps, but the one set by default doesn't suit me at all. Example :

Reference point : default is up and left i need it down and right
Segments are set to 2 left and 4 right, I need 2 right only
Style is simplebox and i need doublebox
...

Do you know a way to default this ? It's a bit tedious to change it manually every time...

For additional info : I didn't find anything in Composer - Settings menu/Composer where i thought i might find those info.


texts: How to improve speed of rendering polygons from Postgis in Leaflet?

texts:
I want to improve speed of rendering polygons in my Leaflet map.

Polygons are stored in Postgis, and are rendered with select query with bbox. I wonder can I cache somehow those polygons?
I read about Tilestache a bit, but I'm still not sure if that would help in speed performance.

Also, one of my problems is that I need to relate postgis tables with non-spatial tables in my postgresql. So, some of the attributes I need to show in the popup are in some other non-spatial tables in Postgresql.

Just to be more clear, that map is already developed and all works but I want to improve speed of polygons rendering. Can you give me some advices?


texts:
I am change the spatial reference system from raster file with gdalwarp. I would like to overwrite file to avoid own the same file again in my computer. I am trying to use gdalwarp with -overite parameter. But the error happened and I didn't understand why. I think that the path of input and output can be the same in order to do what I would like (not have the same file recorded). Is it possible to do with gdalwarp?

gdalwarp -overwrite /tmp/LC08_L1TP_215068_20160913_20170321_01_T1/LC08_L1TP_215068_20160913_20170321_01_T1_B8.TIF /tmp/LC08_L1TP_215068_20160913_20170321_01_T1/LC08_L1TP_215068_20160913_20170321_01_T1_B8.TIF -s_srs EPSG:32624 -t_srs EPSG:4674
ERROR 5: Source and destination datasets must be different.


texts: Problem with gdal in C# to reproject raster

texts:
I wrote this simple code below in C# to convert a raster from Lambert Azimutal Equal Area to regular long-lat raster using GDAL. It works quite fine excepted for the zero-value bands on each side, in purple in the image. I understand that those area would have the nodata_value (i.e. -9999) in turn of zero.

Using comments to this question I found that using a .vrt file and adding this to the   <GDALWarpOptions> section:

<Option name="INIT_DEST">-9999</Option>

it works!  So my problem now is how to add this option programmatically, intended that I'm using a AutoCreateWarpedVRT.

            Dataset ds = Gdal.Open(file, Access.GA_ReadOnly);
            SpatialReference sss = new OSGeo.OSR.SpatialReference(null);
            sss.ImportFromProj4("+proj=laea +lat_0=52 +lon_0=10 +x_0=4321000 +y_0=3210000 +ellps=GRS80 +units=m");
            string srs_wkt;
            sss.ExportToWkt(out srs_wkt);

            SpatialReference ttt = new OSGeo.OSR.SpatialReference(null);
            ttt.ImportFromProj4("+proj=longlat +datum=WGS84");
            string trs_wkt;
            ttt.ExportToWkt(out trs_wkt);

            Dataset tgs = Gdal.AutoCreateWarpedVRT(ds, srs_wkt, trs_wkt, ResampleAlg.GRA_NearestNeighbour, 0);

            Driver ddf = Gdal.GetDriverByName("AAIGrid");
            ddf.CreateCopy(@"out.asc", tgs, 1, null, null, null);


texts: How to restrict features that users can see by area?

texts:
Basically I have around 10 tables that are worked by several teams of users. What I need to avoid is team A complaining that team B got a better area to work on or worse, having them sabotaging each other
To delineate these areas, I have a polygon layer. Ideally, users would only see features in all tables that fall inside this polygon's area and anything outside would not be shown to them
What is the best way of achieving this? layer level security? bounding boxes tied to logins? creating views for each table for each team?


texts: Clipping a raster file into the multiple shapefiles

texts:
I split the single shapefile into the multiple shapefile based upon their attributes. Now I want to clip the raster file into the each split shapefile. How can I do it at once?

I am using ArcGIS version 10. I used spilt layers by attributes tools to split the shapefiles.


texts: Use Python binding of osgeo - Union two linestring layers into one

texts:
I am using the Python GDAL/OGR bindings.  I have two linestring shapefile layers.  They are the census TIGER streets and a new set of line segments that represent lines connecting census block centroids to those streets.  I want to merge the two into one shapefile (so I can do some route analysis) where the connector line segments break the street map segments where they intersect.  I can do this with QGIS but need to automate this over thousands of datasets.

EDIT:

I believe I have GDAL 2.1 on a MacOS (Sierra)

  $ geos-config --version
  3.5.0
  $ gdal-config --version
  2.1.0
  $ gdal-config --dep-libs
  -L/Users/me/anaconda3/lib -lgeos_c -lsqlite3 -lexpat -L/Users/me/anaconda3/lib -lxerces-c -lpthread -L/Users/me/anaconda3/lib -lhdf5 -L/Users/me/anaconda3/lib -lmfhdf -ldf -L/Users/me/anaconda3/lib -lgeotiff -L/Users/me/anaconda3/lib -ltiff -L/usr/lib -lpq -lz -lpthread -ldl -L/Users/me/anaconda3/lib -lcurl -lssl -lcrypto -lssl -lcrypto -lz -liconv -lxml2 -lz -lpthread -licucore -lm -L/Users/me/anaconda3/lib -lkea -L/Users/me/anaconda3/lib -lhdf5 -lhdf5_hl -lhdf5_cpp

The below results in no errors but with an empty Shapefile with 0 features.

Edit 2:  I do get an error.  Segmentation fault: 11 on the union operation which is not an uncommon problem.  I have verified the shapefile created using my CreateShapeFromCSV in QGIS and all seems fine.

def CreateNewShapeLayer(shapeName):
  driver = ogr.GetDriverByName("ESRI Shapefile")
  data_source = driver.CreateDataSource(shapeName)

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

  # create the layer
  layer = data_source.CreateLayer(shapeName.rsplit("/", 1)[-1], srs, ogr.wkbLineString)
  print ("Created layer {}".format(shapeName))

  return layer

def UnionBlockCentroidStreetLines():

  censusstreetlayersrc = "/Users/cthomas/Development/Data/spatial/Network/streets/tl_2016_06000_roads_la_clipped.shp"
  censusStreets = ogr.Open(censusstreetlayersrc, 0)
  censuslayer = censusStreets.GetLayer()

  connectorlayersrc = "/Users/cthomas/Development/Data/spatial/Network/streets/street_segment_block_centroid_connectors.csv"

  # Convert my CSV connector layer to a Shapefile
  CreateShapeFromCSV(connectorlayersrc)

  connectorStreets = ogr.Open(connectorlayersrc.replace(".csv", ".shp"), 0)
  connectorLayer = connectorStreets.GetLayer()

  # Create a new empty shapefile to receive the results of my union
  outputLayer = CreateNewShapeLayer(censusstreetlayersrc.rsplit("/", 1)[0] + "/la_streets_with_block_centroid_connectors.shp")

  print("About to union the layers")

  censuslayer.Union(connectorLayer, outputLayer)

  outputLayer = None
  censuslayer = None
  connectorLayer = None

UnionBlockCentroidStreetLines()

I see I can do something feature by feature, but was hoping I could do it for the whole layer.

The two layers are represented in the image below.


texts: ArcGIS Spatial Join - one-to-many, polygon-to-polygon with simplified boundaries

texts:
DATA:
I have two polygon datasets that have no attributes in common other than location, so must be joined using a spatial join. One dataset is polygons of cadastral parcels (i.e. properties) and the other of hunting districts. Every part of the country is officially assigned to a hunting district, so both datasets completely cover the same area.

DESIRED OUTPUT:
Each cadastral parcel should have all of the hunting districts responsible for it listed. This relationship is clearly one-to-many (i.e. one cadastral parcel to many hunting districts) -- most parcels belong to only one hunting district, but larger parcels often belong to 2+.

PROBLEM:
The problem is that the boundaries between the datasets don't match completely (i.e. the hunting dataset boundaries were slightly simplified at some point) and often overlap most neighboring cadastral parcels ever so slightly, rendering a spatial join based on "INTERSECT" useless. The "HAVE CENTER IN" works perfectly for all parcels that belong to only one hunting area, which is the majority, but produces false results for parcels belonging to more than one hunting area.

This whole analysis would be easier to do in any number of other spatial databases with SQL, but I am unfortunately limited to ArcGIS for this task.

QUESTIONS:
The only thing I can think of is to try a negative search radius.
One can set a negative search radius for "INTERSECT" and "WITHIN A DISTANCE", but I am uncertain as to how these two queries differ in ArcGIS when giving a negative search distance. Do they not then become synonymous?


texts: Loading table in ArcMap gives Serious Application Error?

texts:
I am an ArcGis 10. ( ArcEditor licence) user and I have a problem when trying to load data with the Load Objects Tool in ArcMap.

I want to load a table, but after I choose the input data (equipment.dbf) - Add- and press next , ArcMap instantly shuts down and the well known window appears:

  ArcGis Desktop has encountered a serious application error and is
  unable to continue.

One more thing, I am trying to load the table into a file geodatabase with a geometric network and relationship classes (but that shouldn't be a problem as support says).

Does anyone know why, and what is the solution?


texts: QGIS 2.18.1 - Batch adding raster layers with Python

texts:
I create action, which adding raster from file (for a grid vector layer).
Now, I'll want to use python for batch adding raster layer for selected items.
How can I do it?


texts: Polygon intersection using QGIS virtual layers

texts:
I'm trying to intersect polygon geometries in QGIS by using a virtual layer:

SELECT
    sbqry.rowid AS gid,
    sbqry.geom
FROM
    (SELECT
        ST_Intersection(land_parcels.geometry, flood_zone.geometry) AS geom
    FROM land_parcels, flood_zone
    WHERE ST_Intersects(land_parcels.geometry, flood_zone.geometry) AND NOT ST_Touches(land_parcels.geometry, flood_zone.geometry))
    AS sbqry;

Unfortunately, sbqry.rowid AS gid returns NULL instead of auto-incrementing values.

Does anybody know how to create an unique identifier column named 'gid'? As far as I know virtual layers are based on SQLite/Spatialite.


texts:
I'm currently stuck in my FME Workspace script. Ideally, I'm reading from a file geodatabase and will update the SDE database in SQL server. I have CRC calculators to check for the changes in the Attributes and Coordinates and if the 2 CRC values are different then i will use the data from the File GDB to update the  SDE database.

This is a screenshot of my workspace:

The issue that i am encoutering is it always updates the data because the geometries being stored in the SDE database is not the same with the geometries in the File GDB.

I modified the SDE writer feature type to match the X,Y Origin and scale  from the File GDB feature class and disabled the "Simplify geometry" but the coordinates stored in the SDE Database are still changing.

These are the coordinates from the File GDB and the coordinates stored in the SDE Database:

Also both the SDE and File GDB have the same spatial reference.

I am quite stumped about this. Tried alot of things but it still automatically changes the coordinates when saved to the SDE database.

Any ideas will absolutely be helpful.


texts:
In ArcGIS Desktop we are able to create a file geodatabase and also a feature class.

Is this possible in QGIS?


texts: Trying to clip many layers using Python in ArcMap, getting error that an output already exists

texts:
I have a nationwide shapefile I want clipped into regions. There's not so many regions that I need to script iterating through each region, but there are enough where I want to batch process this. I explored using the Python interpreter built into ArcMap 10.5.1 to do so.

First I tested the process with one region and the following code works:

arcpy.Clip_analysis(in_features="WBDHU12", clip_features="07", out_feature_class="F:/GIS/07_HUC12s.shp", cluster_tolerance="")

It returns:

<Result 'F:\\GIS\\09_HUC12s.shp'>

Then, following instructions on entering multiple commands using the ArcMap Python window, I tried to do two clips for two different regions with one script execution:

arcpy.Clip_analysis(in_features="WBDHU12", clip_features="08", out_feature_class="F:/GIS/08_HUC12s.shp", cluster_tolerance="")
... arcpy.Clip_analysis(in_features="WBDHU12", clip_features="09", out_feature_class="F:/GIS/09_HUC12s.shp", cluster_tolerance="")
...

I'd expect that just as I could clip HUC12s to region 7 using that code successfully, I could clip HUC12s to region 8, and then again to region 9 one after the other. If this worked, I'd just setup the script to run all regions of interest and then let it geoprocess them all at once ('set it and forget it', I hoped). However, I hit this error:

Runtime error  Traceback (most recent call last):   File "<string>", line 1, in <module>   File "c:\program files (x86)\arcgis\desktop10.5\arcpy\arcpy\analysis.py", line 62, in Clip     raise e ExecuteError: ERROR 000725: Output Feature Class: Dataset F:\GIS\08_HUC12s.shp already exists.

I've tried a few different regions like this and anytime I run one region at a time, it takes a couple of minutes then completes the clip, but if I run two regions at a time, it instantly gives me this error and oddly enough it also instantly successfully clips the first region (e.g. in the example above, I'll have a new file in my catalog and map for HUC12s clipped to region 8). What I don't understand is:

Why is the clipping so fast when trying to do 2 regions in one script execution, whereas when doing 1 region per execution it takes longer?
Why am I getting an error about the out_feature_class already existing when I have set different out_feature_classes for each region / instance of clipping?


texts: Reading ArcSDE data of a specific version in Python

texts:
How can I read a version within Python?

For example, let's say I set up a workspace which points to a .SDE file connection. Now how can I read the features from a version within Python? I have looked online but can't find anything. Ideally if there is a code that looks for the version that is currently active and in use within the mxd?


texts: Using multipolygon or polygon features for large intersect operations?

texts:
I have two large layers that I want to do intersect operations. One is municipalities and the other is landuse. Right now, the landuse layer is of MultiPolygon type that is, for the entire state, I have only one entry in the table with the class "Forest", but composed of several small polygons.

If I wanted to find all the forest areas in my city, what would be more efficient in PostGIS:

Do an intersect between my municipality and the FOREST multipolygon that takes up the entire state or,
Convert the Forest multipolygon into a new layer with several polygons and then do the Intersect

My gut feeling is that if I break the multipolygon feature into single polygons, the intersect might be faster because it would only occur in the polygons that have a bounding box that overlaps my municipality. If I have multipolygon features, the intersect would have to go through the entire state to identify only the parts that fall inside the municipality.

This is more or less what I'm trying yo get. I want the light green polygons.

Does it matter it the forest polygons (F) are stored as multipolygons or singlepolygons?

I'm going to be processing a very large dataset (Amazon)


texts: Adding row with default values to new created table with ArcPy?

texts:
In my code I am creating a table with some fields.
The fields should be filled with 0 (Datatype of the fields are LONG or DOUBLE). I want to fill them because afterwards I want to update the values. But without a row with default values, I can't update the values.

    tb_path = r'C:\Users\a48384\Desktop\00_ArcGis\Daten'
    tb_name = '\netzplan'
    arcpy.CreateTable_management(tb_path,tb_name)
    tbNetzplan = r'c:\users\a48384\desktop\00_arcgis\daten\netzplan'

    #names of the fields which will be added
    fn_LEN_NS_FREI = "LEN_NS_FREI"
    fn_NS_VERKAB_GRAD = "NS_VERKAB_GRAD"
    fn_LEN_NS_KAB = "LEN_NS_KAB"
    fn_ANZ_NS_KVS = "ANZ_NS_KVS"
    fn_ANZ_NS_KVS_BS = "ANZ_NS_KVS_BS"
    fn_ANZ_NS_KVS_NBS = "ANZ_NS_KVS_NBS"
    fn_ANZ_NS_KVS_NULL = "ANZ_NS_KVS_NULL"

    #Adding fields to the table
    arcpy.AddField_management(tbNetzplan, fn_LEN_NS_FREI, "Double")
    arcpy.AddField_management(tbNetzplan, fn_NS_VERKAB_GRAD, "INTEGER")
    arcpy.AddField_management(tbNetzplan, fn_LEN_NS_KAB, "Double")
    arcpy.AddField_management(tbNetzplan, fn_ANZ_NS_KVS, "INTEGER")
    arcpy.AddField_management(tbNetzplan, fn_ANZ_NS_KVS_BS, "INTEGER")
    arcpy.AddField_management(tbNetzplan, fn_ANZ_NS_KVS_NBS, "INTEGER")
    arcpy.AddField_management(tbNetzplan, fn_ANZ_NS_KVS_NULL, "INTEGER")

Now I would like to add a row with default values. But I don't know how... Afterwards the values are going to be updated:

    cur = arcpy.UpdateCursor(tbNetzplan)
    for row in cur:

        row.insertValue(fn_LEN_NS_FREI, VALUE1)
        row.setValue(fn_NS_VERKAB_GRAD, VALUE2))
        row.setValue(fn_LEN_NS_KAB, VALUE3)
        row.setValue(fn_ANZ_NS_KVS, VALUE4)
        row.setValue(fn_ANZ_NS_KVS_BS, VALUE5)
        row.setValue(fn_ANZ_NS_KVS_NBS, VALUE6)
        row.setValue(fn_ANZ_NS_KVS_NULL, VALUE7)
        cur.updateRow(row)


texts: I'm building a webGIS using leaflet library. Since the data visualization is for only a certain region, it seems that the broader leaflet base map looks a bit distracting. I wonder if there's a way to hide the redundant region.


texts: How to find buffer radius around raster cell that covers a defined sum of raster cells?

texts:
I want to do a neigbourhood analysis on some raster data. I need to know how big a circle around each raster cell must be to include a certain amount of specified raster cells.

Let's say that for both rasters I want to know the radius of the circle needed to reach an aggregated value of 5. In the first case the "circle" must have a radius of 2 to reach the sum of the desired 5 because a radius of 1 would only include one raster cell that I am interested in. In the second case already a smaller circle with a radius of 1 is sufficient.
My outcome should thus be a raster file where every cell has the value of the radius that is needed to reach that defined sum.
I was thinking a long time about a good title for this question. I would welcome better suggestions.


texts: Is there an open source stereoscopic interface for QGIS, like PurVIEW or DAT/EM, or stereo Analyst?

texts: I have a lot of imagery processing that I do at the local college, and most of the imagery is stereo imagery.  Is there an open source piece of software that could be used to view the imagery and perhaps integrate it with open source GIS software, like QGIS?


texts: Saving bookmarks in a database using javascript

texts:
I want to save a bookmark permanently in a database using javascript. I have an idea.

"I have a service. so that i shall use that service to interact my code and the database."

I am struck with this how to begin. I would like you to suggest whether this is possible. If so how to proceed with this?


texts:
I have done point layer overlay with OpenLayers. Whenever I click any point, it's showing me all the data of the table (PostgreSQL) in popup.

How to control amount of data in the popup?

All the columns are being displayed in popup. I just want location_name to be displayed in the popup. How to extract location_name from event.text?

Code is:

....
info = new OpenLayers.Control.WMSGetFeatureInfo({
    url: 'http://localhost:8080/geoserver/layername/wms',  // geoserver point layer
    title: 'Identify features by clicking',
    queryVisible: true,
    eventListeners: {
        getfeatureinfo: function(event) {
            map.addPopup(new OpenLayers.Popup.FramedCloud(
                "ABCD",
                map.getLonLatFromPixel(event.xy),
                null,
                event.text,
                null,
                true
            ));
        }
    }
});
map.addControl(info);
info.activate();


texts: OpenLayers with a small, discrete set of XYZ tiles

texts:
I'm trying to set up a handy map of my Minecraft server to share with my fellow players; the specifics on that front aren't really necessary though -- suffice it to say that I have a tile set (each tile 256x256 pixels) representing 6 zoom levels, with a 2x2 grid of tiles at zoom level 0 all the way up to 64x40 tiles at zoom level 5.

After much tinkering and excessive frustration, I have managed to finally get the map fully viewable -- but OpenLayers thinks I've got tiles that are outside my bounds, leading to the map being off-center when initially loaded and, more importantly, leaving me thinking that I've merely gotten lucky with the OpenLayers.Bounds declaration and that it will all fall apart completely when the map (and its dimensions) are next updated.

Here's what I'm doing to instantiate the map (additional commented-out OpenLayers.Bounds declarations show some of my thought processes and frustrations):

var map; //The map
var day, night; //The layers, one for a "day" map and one for a "night" map

function init(){
    map = new OpenLayers.Map('map', {
        tileSize: new OpenLayers.Size(256,256), //All my tiles are 256x256
        //I've tried to many different bounds...
        maxExtent: new OpenLayers.Bounds(0,0,16384,10240), //Size of my map at zoom level 5
        //maxExtent: new OpenLayers.Bounds(0,0,512,512),
        //maxExtent: new OpenLayers.Bounds(0,0,8138,5034), //Size of the raw .png image
        //maxExtent: new OpenLayers.Bounds(0,0,8192,8192), //Size of the raw .png image, rounded up to a multiple of 256
        //maxExtent: new OpenLayers.Bounds(0,0,3000,2500),
        maxResolution: 16, //Anything else seems to be even worse -- what's so magical about 16?
        numZoomLevels:6, //Cuz that's how many tile sets mcmap generates
        units: "pixels", //Not listed in the docs, but found on another site -- legit?
        projection: "", //Not listed in the docs, but found on another site -- legit?
        eventListeners: {
            "zoomend": mapEvent,
            "moveend": mapEvent
        }
    });

    //Add the "day" layer
    day = new OpenLayers.Layer.XYZ("Day", "./tiles/x${x}y${y}z${z}.png");
    map.addLayer(day);
    //"Night" layer omitted for brevity

    map.zoomToMaxExtent();
}

This displays my 2x2 grid of tiles at zoom level 0 fine, but in the corner of a larger 4x3 grid -- and I can't figure out why! (As of the time of this writing, a live demo of what I'm getting can be found here.)

I believe my problem is with the maxExtent and/or maxResolution options; frankly, I can't make heads or tails of how either of these values actually relate to the tiles I have, especially since "16" seems to be the magic number for the latter, for no reason I can see!

How do I go from a simple set of tiles that I want arranged on an XYZ grid, to a functional OpenLayers map, without OpenLayers requesting extraneous bogus tiles outside the confines of my tile set? Merely having the right values for this particular tile set is insufficient -- my tile set will change over time as more of my server's map is explored and revealed, so I need to be able to determine the appropriate values/settings for any arbitrary tile set (well, with 6 zoom levels, at least).

Edit: Just to clarify, I know how to make the pink tiles "disappear", that's not the question -- the question is how do I convince OpenLayers that those tiles are "out of bounds" and not even request them in the first place?


texts: Lines on reprojected sp objects with Mollweide projection

texts:
I'm having an odd problem such that when I reproject maps to a Mollweide projection in R using sp, I get additional lines stretching across the globe. What is going on here, and do you have any suggestions to fix it.  Here's a reproductible example:

library(sp)
library(maps)
library(maptools)

worldmap <- map("world", plot=F)
plot(worldmap, type="l)

So far, so good. I get

So now I reproject. This is a CRS I'm using to be compatible with some other projects. Note how world-spanning lines then crop up.

worldmapLines <- map2SpatialLines(worldmap,
                                  proj4string=CRS("+proj=longlat +datum=WGS84"))

worldM <- spTransform(worldmapLines,
                      CRS("+proj=moll +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +units=m +no_defs"))

plot(worldM)


texts:
My ArcPy script doesn't run because a syntax error was raised, but I already had checked the syntax and indentation of the entire script.

Can anyone help me?

The error that was raised points to a colon that's the end of 'if' on the line 236 as you see in the print screen below.


texts: Label callout lines to the closest point on the label

texts:
I am trying to improve callout lines on maps I'm making. I'm currently using the method in this link to make my callout lines. That works fine, bu I'd like the lines to always go to the closest point of the label like this:

I think it would work if I used code down the line of this in the expression editor for the geometry generator.

make_line(
  make_point($x, $y),
    make_point((CASE
      WHEN "auxiliary_storage_labeling_positionx" > $x
        THEN "auxiliary_storage_labeling_positionx"
      WHEN ("auxiliary_storage_labeling_positionx"  + "labelwidth") > $x
            AND ("auxiliary_storage_labeling_positionx") < $x
        THEN ("auxiliary_storage_labeling_positionx" + "labelwidth")
      ELSE $x
      END )
    ,(CASE
      WHEN "auxiliary_storage_labeling_positiony" < $y
        THEN "auxiliary_storage_labeling_positiony"
      WHEN ("auxiliary_storage_labeling_positiony"  + "labelHeight") < $y
            AND ("auxiliary_storage_labeling_positiony") > $y
        THEN ("auxiliary_storage_labeling_positiony"  + "labelHeight")
      ELSE $y
      END )
  )
)

I'm not sure if I have all the < and > correct but that's beside the point.
The issue with that is that I can't seem to figure out how to get the width and the height of the label.

I am not very good with the expression editor, and have never done anything with the function editor.


texts: How to connect Qgis with pycharm in Ubuntu OS

texts:
I have just started using pycharm with qgis but could not connect both of them. Pycharm always stays in a state of "waiting for connection". Most the tutorials available points to the windows but I am using ubuntu so could not find way to debug qgis code on pycharm. Here is my pycharm code:

from shapely.geometry import *
from shapely.wkt import loads

import sys

import pydevd

pydevd.settrace('localhost', port=53100, stdoutToServer=True, stderrToServer=True)

class Loader:

    def __init__(self, iface):

        """Initialize using the qgis.utils.iface
        object passed from the console.

        """
        self.iface = iface

I have enabled breakpoints in the pycharm and also added the pycharm-debug.egg in the pythonpath does anybody how configure it from qgis on ubuntu?

pycharm is always in:

Starting debug server at port 53100
Use the following code to connect to the debugger:
import pydevd
pydevd.settrace('localhost', port=53100, stdoutToServer=True, stderrToServer=True)
Waiting for process connection...

When I run this script on the top from qgis nothing happens, the breakpoint does not get called.


texts:
I'm trying to calculate the direction the 'generate near table' function gives me into actual directions, N, NW etc.  Azimuth would be a start too.
I tried this script in the field calculator but it doesn't run.

This was the original we tried:

if([Direction]>0,
    if([Direction]<22.5,"W",
    if([Direction]<67.5,"SW",
    if([Direction]<112.5,"S",
    if([Direction]<157.5,"SE",
    if([Direction]<180,"E",""))))),

if([Direction]>-22.5,"W",
    if([Direction]>-67.5,"NW",
    if([Direction]>-112.5,"N",
    if([Direction]>-157.5,"NE",
    if([Direction]>-180,"E",""))))))

I tried it with and without "then output=" and "Else, output =0".
I'm using ArcGIS 10, have tried with with python and vb ticked.


texts: How to create tiles (compatible with OSM) from a MapServer mapfile

texts:
For an internship, I have to make an S57 online map server which will be used later for something else.
I have some constraints:

My cards must be OSM compatible (the existing version of the final product using an OSM card)
I must be able to have geographical references, be able to display a point of given latitude and longitude.

From a previous project I have a pretty good understanding of the S57 standard.
So I started looking for a utility that could display S57 charts online. I came across MapServer (with .map files) and this code (https://github.com/LarsSchy/SMAC-M) first. The SMAC-M code generates .shp and .map files which are used by MapServer to display maps.
So already I have the ability to display the maps in a browser.
But the result is not compatible with OSM; to have one that suits me, I chose to use MapCache to convert MapServer's output to OSM equivalent tile.
Here is the discussion I had on github
https://github.com/LarsSchy/SMAC-M/issues/27
I want to know if is the best choice?
If yes:
I would look for sample configuration code for MapCache that does what I want to do.
And also in case of tutorials on MapCache (there is 1 in the documentation but it is not sufficient).
Finally, I want anyone can help me write an ideal configuration for what I want to do.
I have tried to modify the MapCache quickstart XML file.
the original file was :
<?xml version="1.0" encoding="UTF-8"?>

<mapcache>
  <cache name="disk" type="disk">
    <base>/var/cache/mapcache</base>
    <symlink_blank/>
  </cache>

  <default_format>PNG</default_format>

  <source name="lake_source" type="wms">
    <getmap>
      <params>
        <FORMAT>image/png</FORMAT>
        <LAYERS>lakespy2</LAYERS>
        <MAP>/usr/local/www/docs_maps/mapserver_demos/itasca/itasca.map</MAP>
      </params>
    </getmap>
    <http>
      <url>http://localhost/cgi-bin/mapserv?</url>
    </http>
  </source>

  <tileset name="lakes">
    <source>lake_source</source>
    <cache>disk</cache>
    <grid>GoogleMapsCompatible</grid>
    <format>PNG</format>
  </tileset>

  <service type="wms" enabled="true">
    <full_wms>assemble</full_wms>
    <resample_mode>bilinear</resample_mode>
    <format>PNG</format>
    <maxsize>4096</maxsize>
    <forwarding_rule name="catch all">
      <http>
        <url>http://localhost/cgi-bin/mapserv?map=/usr/local/www/docs_maps/mapserver_demos/itasca/itasca.map</url>
      </http>
    </forwarding_rule>
  </service>
  <service type="wmts" enabled="true"/>
  <service type="tms" enabled="true"/>
  <service type="kml" enabled="true"/>
  <service type="gmaps" enabled="true"/>
  <service type="ve" enabled="true"/>
  <service type="mapguide" enabled="true"/>
  <service type="demo" enabled="true"/>

  <log_level>debug</log_level>

  <errors>report</errors>
  <locker type="disk">
    <directory>/tmp</directory>
    <timeout>300</timeout>
  </locker>

</mapcache>

I have change mapfile path with my own mapfile path.
but the result was the same.


texts: OpenStreetMap- Field Calculator to Merge/Flatten Layer Attributes?

texts:
I'm not sure exactly what I am asking, probably lending to my struggle to figure it out via google.

I have a large polygon layer from OpenStreetMap data (openstreetmap.org). My problem is I need to reduce/generalize the fields/attributes in that layer largely due to inconsistencies in classifications from different people adding to it.

For example, I have the value field "leisure" with polygons labeled "garden", "golf course", etc. I have another value field with polygons/attributes "grass", "meadow", etc.

Is there a way to create a new value field or separate polygon layer with all these polygons labeled as "greenspace"?

I am still struggling to figure out how to use the field calculator to do what I am looking for. I think my hang up is that I am trying to combine polygons of different "attributes" within different field types (attaching image to original post for clarification). Is there a way to combine polygons of several values (not fields) from different fields into one field?


texts: Map can only have one feature, after user click

texts:
I have a map with a vector layer and Editing tools. I want the map to have only one feature even if the user keeps adding. So I want the new feature to automatically erase the previous one (if a previous exists). That way there is only one feature in the map, the new one that the user adds (using Editing tools)

I tried some things, but I end up erasing the feature, as soon it is added. How can I do this?

Here is my code:

myLayer.events.register('featureadded', this, feature_added);
map.addControl(new OpenLayers.Control.EditingToolbar(myLayer));

function feature_added(feature){
   //erase previous features here (?)
   pins.events.unregister('featureadded', this, feature_added);
}

I'm using OpenLayers 2.13.1.


texts: QGIS.2.18. cannot save calculated area into layer

texts:
After I calculated several areas in the attribute table (added a column, named it "Areas" and then calculated with "$area" in the calculator), I could not save the new information within my given layer.

Instead I keep getting the below message
"Could not commit changes to layer landuse
"Errors: ERROR: 1 new attribute(s) not added
  ERROR: 55 attribute value change(s) not applied.
  Provider errors:
      SQLite error: duplicate column name: AREA
    SQL: ALTER TABLE "landuse" ADD COLUMN "AREA" FLOAT
      SQLite error: near "WHERE": syntax error SQL:"

How to deal with that?

I am very new to QGIS and have no coding background.


texts: Is there a way to identify and remove similar GPS tracks in PostGIS or ArcGIS?

texts:
I have a large number of GPS tracks from a number of different users. I'd like to identify tracks that take the same route and remove those duplicates. For example, a user takes the same streets from the same origin/destination, but obviously with GPS accuracy issues, the lat/longs aren't exactly the same.

Any sort of intersect function simply provides a large set of intersecting sections, but doesn't compare a full track from start to finish.

Currently using PostGIS, but trialing in ArcGIS to visualize and test a procedure before applying it to the 1,000+ tracks and 8 million datapoints.


texts: Failed to create address locator from a CSV file

texts: I am using ArcMap Desktop 10.4. I have a CSV file containing the address of properties in 2 fields - civic number and street name. I was trying to create an address locator using the file as the reference data. I selected single house subaddress and made sure all the fields match. But when I clicked OK at the end, it just says failed to create the address locator(Error 000042). I've cleaned the file and tried a few times but it was the same result. what are some possible explanations and how can i work around that? Without creating the address locator, I can't geocode.


texts: Which is the correct Options 's "Hydraulics" field in EPANET?

texts:
Which parameter has to be given into Options 's "Hydraulics" field in EPANET?
Has it to be only a name or a complete file name with the extension?
In my case .rpt doesn't work, neither without the extension (nor with a path).
[I am working with giswater, that generate new db in postgreSQL, files for EPANET and a .qgs]


texts: QGIS export to GPX: Error "Geometry type of 'multi polygon' not supported" (GPS tools)

texts:
I am trying to convert a shapefile with multiple polygons to a GPX file in QGIS 3.10.3, which has now repeatedly failed with varying export settings since "MultiPolygon" is not supported.

I have already merged all the different polygons in the shapefile, so I only have one polygon in the attribute table. Still I get this error code. Any tips?

Thanks in advance!

Elena

PS: my final goal is to import these shapefiles to my Garmin Etrex Vista HCx.


texts:
I want to work with landsat 7 images for land cover and NDVI evaluation. There are strips in the landsat images, and I don't know how to remote them. I have learned methods for gapfill but I don't know how to work with. I also used Erdas 2011 (focal analysis) for the processing but I see no result. I have also a sort of gap mask download with my landsat 7 images but I don't know how to use it.

How should I deal with the striping?


texts: Is there a way to get the X,Y from a feature class?

texts:
I currently have an instance of a selected feature class, and I want to be able to get the midpoint X,Y (if its polygon) or the X,Y (if its IPoint) of the feature from the layer, is there some build in GIS method that I can achieve that easily?

Thanks and Regards,


texts: Setting the number of bars in a histogramm

texts:
I am using qgis' 2 "Statist" to build a histogramm for my data. I see that I can set the band width of the bars, but it seems that "Statist" does not provide the option of setting the number of bars.
Any ideas on how to set the number of bars in a histogramm would be much appreciated!


texts: Unclear on proxy used for AJAX in Javascript

texts:
I do understand how the proxy helps in transferring the XMLHttpRequest between the ArcServer and the Server (IIS/Apache) but where I am confused is that the first time the Html Doc loads, doesn't the Server need to request the ArcServer then as well. How does that request complete without the proxy?

I mean the first time the page loads it does not need that proxy and still the map is loaded in the Map container. But later if we send a XMLHttpRequest to the Geoprocessing Tool, the response from the ArcServer can only be received by the Server if we have setup the proxy, right?

So how does it happen the first time and not the second time?


texts:
Is there a way to set up an automated photo import from a network folder in QGIS?

Specifically, I would like to push a custom button and select a root folder, and as images populate the folder I would like them to import into QGIS.  Intervals of 30 min to 1 hour for automated import but with the option of manually refreshing and importing if there are new files.

Is this even possible? Not looking for hand holding or completed answer, maybe a nudge in the right direction.


texts: font error when using mapnik generate_tiles.py

texts:
generate_tiles.py has worked fine throughout the transition from mapnik 0.7 to mapnik 2.x in gentoo linux.  Now I have switched to debian (wheezy) and I am getting the following error when I run generate_tiles.py:

RuntimeError: Failed to find font face 'unifont Medium' in FontSet 'book-fonts' in map '.../osm.xml'

If I comment out the book-fonts section in fontset-settings.xml.inc, I get:

RuntimeError: Failed to find font face 'unifont Medium' in FontSet 'bold-fonts' in map '.../osm.xml'

However, using tilelite creates the tiles just fine.

$ ls `python -c "import mapnik2;print mapnik2.fontscollectionpath"`
total 9048
drwxr-xr-x  2 root root   4096 2014-10-11 09:10 .
drwxr-xr-x 13 root root   4096 2014-10-11 20:48 ..
-rw-r--r--  1 root root 611212 2011-02-27 07:37 DejaVuSans-BoldOblique.ttf
-rw-r--r--  1 root root 672300 2011-02-27 07:37 DejaVuSans-Bold.ttf
-rw-r--r--  1 root root 580168 2011-02-27 07:37 DejaVuSansCondensed-BoldOblique.ttf
-rw-r--r--  1 root root 631992 2011-02-27 07:37 DejaVuSansCondensed-Bold.ttf
-rw-r--r--  1 root root 576004 2011-02-27 07:37 DejaVuSansCondensed-Oblique.ttf
-rw-r--r--  1 root root 643852 2011-02-27 07:37 DejaVuSansCondensed.ttf
-rw-r--r--  1 root root 345204 2011-02-27 07:37 DejaVuSans-ExtraLight.ttf
-rw-r--r--  1 root root 235848 2011-02-27 07:37 DejaVuSansMono-BoldOblique.ttf
-rw-r--r--  1 root root 313856 2011-02-27 07:37 DejaVuSansMono-Bold.ttf
-rw-r--r--  1 root root 241972 2011-02-27 07:37 DejaVuSansMono-Oblique.ttf
-rw-r--r--  1 root root 333636 2011-02-27 07:37 DejaVuSansMono.ttf
-rw-r--r--  1 root root 611556 2011-02-27 07:37 DejaVuSans-Oblique.ttf
-rw-r--r--  1 root root 720012 2011-02-27 07:37 DejaVuSans.ttf
-rw-r--r--  1 root root 332036 2011-02-27 07:37 DejaVuSerif-BoldItalic.ttf
-rw-r--r--  1 root root 341072 2011-02-27 07:37 DejaVuSerif-Bold.ttf
-rw-r--r--  1 root root 331128 2011-02-27 07:37 DejaVuSerifCondensed-BoldItalic.ttf
-rw-r--r--  1 root root 316440 2011-02-27 07:37 DejaVuSerifCondensed-Bold.ttf
-rw-r--r--  1 root root 338140 2011-02-27 07:37 DejaVuSerifCondensed-Italic.ttf
-rw-r--r--  1 root root 330012 2011-02-27 07:37 DejaVuSerifCondensed.ttf
-rw-r--r--  1 root root 338776 2011-02-27 07:37 DejaVuSerif-Italic.ttf
-rw-r--r--  1 root root 363200 2011-02-27 07:37 DejaVuSerif.ttf

No book font, but there are bold fonts.

These are all the dejavu fonts installed (with their version numbers):

ttf-dejavu                                          2.33-3
ttf-dejavu-core                                     2.33-3
ttf-dejavu-extra                                    2.33-3

I've seen a few people having this problem around the web, but no solutions offered.


texts: CSV table import: All rows from CSV file represented in one row in ArcMap attribute table

texts:
When I import the following csv file, all rows only appear in one line in attribute table with the headings as "Field#".

I copied original excel export in new sheet where I pasted values only. I've never had this problem before with importing csv files in ArcMap. Any ideas?


texts: Rotating Polygon in 3D Space as KML for Google Earth?

texts:
I'm trying to create this effect in a KML to be displayed in Google Earth.  I am writing in KML 2.1 as a reference.

I know it is possible to create a polygon that is circular in appearance, but can that same polygon be "stood on its edge," similar to how a U.S. penny sits when it is showing both the heads and tails side?


texts: Can we create a temporary in-memory geodatabase using python?

texts: I have been using in-memory workspace for creating temporary features, but I was wondering if we could create a temporary in-memory geodatabase (using python) and then use it for importing other in-memory features?


texts:
I'm trying to add a PDF into my ArcGIS. As I understand it is possible to change the PDF to a DWG, but I don't know how to do this.

Any ideas?


texts: Link ui.Thumbnail to external URL in Earth Engine?

texts:
I figured out how to import logos into Earth Engine by following the link below..
Adding a logo to a panel on an app in Google Earth Engine?
I would now like to code the logo such that when it is clicked, it will open a new tab from the link provided. I believe this should be possible because onClick functions can be written for the thumbnail.
I have put together the following code...
// DEFINE MAIN APP PANEL
var main = ui.Panel({style: {width: '310px', padding: '8px'}});

// IMPORT LOGO
var logo = ee.Image("users/SE_FHM/Logos/Tall_Timbers_Logo")

// DEFINE WEB LINK FOR THUMBNAIL
var link = ui.Label('Tall Timbers Research Station', {},'https://www.talltimbers.org');

// CONVERT TO THUMBNAIL, INSERT ON CLICK FUNCTION TO OPEN UP THE WEB LINK
var branding = ui.Thumbnail({image:logo,params:{bands:['b1','b2','b3'],min:0,max:255},style:{width:'270px',height:'123px'}});

// ADD TO THE MAP
main.add(branding);
main.add(link)
ui.root.insert(0, main);

Link to full code
Are there any easy solutions to achieve this?


texts:
I am very new in QGIS.

I imported a DEM in QGIS, I generated contours, then I have added some labels.

I have been able to export those contours in DXF but the labels haven't been exported.

How can I export the contours with labels?


texts:
I'm using ArcGIS and closely following this tutorial http://economics.mit.edu/files/8945 to process data in ArcMap. My problem arises because in my software I can't find the Create Fishnet option in Data Management Tools-->Features Class. Here is a picture of what it looks like

I have searched in the Internet but no one seems to be having the same issue. The coordinate system is the standard WGS 1984 Web Mercator.


texts:
I'm having trouble in QGIS with the symbology. I can classify and categorized my data, but I cannot with the graduated symbology?

I've try different format in Excel, and I can't find the problem here.


texts: How to give a value of 0 in a raster only in a certain area?

texts:
I want to give the value of 0 in this raster but only in the area within the polygon. So the final outcome is the same raster file but with 0 values in that area
How can I do that?


texts: Updating multiple fields with list values using arcpy.da.UpdateCursor

texts:
I have a script that contains a function that fetches XML from a web service, containing a series of metadata values.  The function takes a field from a feature class as criteria, queries the web service, and returns a list of metadata values associated with the submitted criteria.  My intent is to update various fields in that feature class with these values, and stuff them into a dictionary so I don't have to make more web calls than necessary.  I am failing at the point where it is trying to update the fields.  Here is the function.  It is not efficient, but it works... when I get everything else working, I'll rewrite it using elementree :)

import datetime, time, urllib, sys, arcpy

def getMetadata(id):
    iDate = ''
    iCategory = ''
    iType = ''
    svcURL = r'http://www.someurl.com/?query=" + str(id)
    localFile = r'C:\Temp\localfile.xml'
    urllib.urlretrieve(svcURL, localFile)
    readFile = open(localFile, 'r')
    for line in readFile.readlines():
        if line.find("<date_time>") > -1:
            iDatePre = line[line.find("<date_time>")+11:line.find("</date_time>")]
            iDateObject = datetime.datetime.strptime(iDatePre, '%d-%b-%y %H.%M.%S.%f %p')
            iDate = datetime.datetime.strftime(iDateObject, '%Y-%m-%d %H:%M')
        if line.find("<Category>") > -1:
            iCategory = line[line.find("<Category>")+10:line.find("</Category>")]
        if line.find("<Type>") > -1:
            iType = line[line.find("<Type>")+6:line.find("</Type>")]
        metaData = [iDate, iCategory, iType]
    readFile.close()
    return metaData

So the idea is to use da.UpdateCursor to put those values in a feature class, using row[0] as the function's criteria:

dictOfMetadata = {}

fgdb_name = 'TESTGDB'
fgdb_path = r'C:\Temp'
fgdb = fgdb_path + '\\' + fgdb_name + r'.gdb'
arcpy.env.workspace = fgdb
fdsb = 'dataset_b'

if arcpy.Exists(fdsb):
    print('Found {} dataset.'.format(fdsb))
    fcs = arcpy.ListFeatureClasses('*','',feature_dataset=fdsb)
    for cfc in fcs:
        with arcpy.da.UpdateCursor(cfc, ['ID','DTG','CAT','TYPE']) as cursor:
            for row in cursor:
                idSend = row[0]
                if dictOfMetadata.has_key(idSend):
                    print('Found local metadata for {}.'.format(idSend))
                    cursor.updateRow(row[1, 2, 3] + dictOfMetadata[idSend][0,1,2])
                else:
                    cursor.updateRow(row[1, 2, 3] + getMetadata(idSend)[0, 1, 2])
        del cursor

The error is: TypeError: list indices must be integers, not tuple

I've tried a bunch of different ways to do this but this is the closest I've gotten.


texts: One (integer) to Many (array of integers) relation in QGIS

texts:
Is it possible to make a One to Many relation betwen two layers where the id of the "One" layer is an integer but the id of the "Many" is an array of integers in postgis? (integer[])

Doing the normal relation in Project -> Properties -> Relations it is not working...

The integer[] field in QGIS is detected as QString, maybe this is the source of the problem...

I am using QGIS 2.18.13


texts: QGIS Field Sequential Numbering with ordering rule from SW to NE

texts:
How could I use the $Rownum tool to get sequential number for a field to make sure the southernmost and westernmost point is 1 and the following numbers follow this rule set (where South is first and West is second importance). For example, when I use $Rownum I have this output:

But What I want is this output: