texts: Error Quantifying Land Cover Classes by Polygon in PostGIS

texts:
I am trying to ultimately calculate proportion of land cover types by polygon in PostGIS, so starting using st_valuecount to get counts of pixels with different land cover types. The data are pretty large (6" resolution land cover data and about 800,000 polygons, varying in size quite a bit).
The code below works as expected issue when applied to a subset of polygons, but when applied to the entire dataset, throws the following error:
SQL Error [XX000]: ERROR: invalid memory alloc request size 1144274944

Here's the code (adapted from PostGIS in Action, 2nd Ed.):
with
    cte as (
        select bbl,
            st_valuecount(
                st_clip(st_union(p.rast), geom_2263),1, true, ARRAY[0,1,2,3,4,5,6,7]
                ) as pv
            FROM
                base_rasters.landcover6in as p
                inner join
                results_scratch.polys
            on st_intersects(p.rast, geom_2263)
            group by bbl, geom_2263
        )
SELECT bbl,
    (pv).value, sum((pv).count)
from cte
group by bbl, (pv).value
ORDER by bbl, (pv).value;

Of note, bringing the exact same datasets into QGIS 3.2.2 and running the 'Zonal Histogram' tool works without issue.
Any suggestions on what might be going on, or how to troubleshoot? (As much as anything, looking to have reproducible SQL code for doing similar tasks in the future.
Here's the software information (running on Win 7 x64 Professional):
PostGIS Version
POSTGIS="2.4.4 r16526" PGSQL="100" GEOS="3.6.2-CAPI-1.10.2 4d2925d" PROJ="Rel. 4.9.3, 15 August 2016" GDAL="GDAL 2.2.4, released 2018/03/19" LIBXML="2.7.8" LIBJSON="0.12" LIBPROTOBUF="1.2.1" RASTER

PostgreSQL Version
PostgreSQL 10.3, compiled by Visual C++ build 1800, 64-bit

Lastly, I'll note I've previously posted about this on the PostGIS listserv but had not received responses, so posting here. I'll share any useful answers reciprocally.
Edit 1:
As noted in the comments, doing this operation without the union results in too many pixels being counted for some polygons (some of the pixels are of classes that do not even lie within the polygon). Thus, doing this without the union would likely work, but I'm observing issues in the clipping of the raster to the polygons, as follows:
Using the following code produces expected result for clipping to a single polygon:
select st_clip(st_union(p.rast), geom_2263)
 FROM base_rasters.landcover6in as p,
      results_scratch.polys
      where st_intersects(p.rast, geom_2263) and bbl=5011720090
      group by geom_2263;

Here's the result - the yellow outline is the a focal polygon (bbl 5011720090); blue lines are borders of the raster tiles that intersect that polygon; the raster seen is the result of the clip, and the darker green filled polygons are the suite of polygons I'm pulling from.

Here's the result when I remove the st_union:
select st_clip(p.rast, geom_2263)
 FROM base_rasters.landcover6in as p,
      results_scratch.polys
        where st_intersects(p.rast, geom_2263) and bbl=5011720090;

Same color scheme as above; the black outside of the yellow/blue boundaries is 0 for nodata. You can see there are areas that are outside of the focal area (yellow boundaries) but within an associated tile that are maintained in this clip.

Any suggestions on what's going on here?
Edit 2
It appears the issue with st_clip across un-unioned raster tiles that I'm dealing with is a documented issue - https://trac.osgeo.org/postgis/ticket/3457 and documented in this stackexchange post.
Of course insights are welcome if I might be missing an alternative way of doing this in PostGIS.
Edit 3
Looks like there's now a bug fix, pending for the next PostGIS release, and with the ticket is now a way to redefine the function to fix for the interim.


texts:
After a frustrating few days I came to the conclusion that the Slippy Map Wiki only works for raster tiles, not for OSM2VectorTiles data. I can work with the data once I find the correct tile_data BLOB, but actually locating the correct tile via a user lat/long is problematic. The equations given in the wiki above identify raster tiles correctly, but I'm trying to get this working for an offline version of OpenStreetMap rendering/navigation (meaning I can't just use mapbox's JS API to convert it, also I can't find where the conversion happens in the github repo).

Can anyone offer any advice as to how I can convert a lat/long to the correct tile in an OSM2VT SQLite database?

The query, that works for raster tiles, may be of some use:

    public void UserQuery()
    {
        Console.Write("Please enter a latitude value: ");
        double userLat = Convert.ToDouble(Console.ReadLine());
        Console.Write("Please enter a longitude value: ");
        double userLong = Convert.ToDouble(Console.ReadLine());

        _merc.TileLat = Math.Floor((1.0 - Math.Log(Math.Tan(DegreesToRadians(userLat)) + 1.0 / Math.Cos(DegreesToRadians(userLat))) / Math.PI) * Math.Pow(2, ZoomLevel));
        _merc.TileLong = Math.Floor((userLong + 180.0) / 360.0 * Math.Pow(2, ZoomLevel));
    }


texts: I've made a nice roadmap and am trying to add a basemap, but the World Topo Map has town names, so does the USA Topo map and the Nat Geo Topo map. how do I get rid of these?


texts: Viewshed analysis using polygons in ArcGIS Desktop?

texts:
I'd like to do a viewshed analysis of mines in the midwest. What I have are point features indicating the location of the mines. I also have an acreage field, which I could use to expand from the point and create a buffer polygon that approximates the size of the mine. I'd like to do a viewshed analysis on this polygon, rather than the point, since the point is just an arbitrary "center" of the mine.

Unfortunately, Viewshed in ArcGIS only allows me to approximate the viewshed of a point itself. Are there any ways to do polygons instead?


texts: Count number of stream segments in a particular stream order: Bifurcation ratio

texts:
Question regarding stream Morphometry.  In which we calculate number of stream segment of a particular stream order.

Issue is that ArcGIS is giving incorrect no of segments. What happened in ArcGIS 10.1 after converting into Vector format (feature class), Grid code shows stream order. But  the no. of segments is different. I mean In this figure we have three segments of  2 order streams [Blue streams] but ArcGIS counted it 7 segments.

Higher order streams got disconnected when new lower order stream join ( at nodes). Although on clicking it shows same order but its  frequency is different i.e  no of segment will be different.

But if two same order join and they meet higher order then stream segment should not be joined ( such as  when two 2nd order ( Blue ones) join with 3rd (higher order) (red ones) in this these two blue streams will not be joined. Because they are different tributary.

If we do this correctly then we can get other Morphometry parameters such as Bifurcation ratio.

Could you help me to sort it out?  We can count manually but for larger area it is not a good way.

I am new to Python coding.


texts:
I have a shapefile of points and want to measure the distance between each of them. After running the Distance Matrix command I get a matrix with distances as expected, but have no idea how to match each row (i.e. each point) to each other.

The problem arises because the names of the columns is not informative of which point each distance refers to.

Clearly the first column refers to the distance between each point and itself. But, for example, how do I know what two points define the distance in row 1, column 2?


texts: Error executing CQL on layer other then demo layers

texts:
I have some sort of wierd problem as I was playing with CQL. CQL works well with given layer but when i try with layer loaded from mapinfo or some views then it does not work. As in the picture when I click the feature it shows the column called NAME  but when I query using  it (or other column) generates error like:

Property &apos;NAME&apos; could not be found in DISTRICT_BND

Why is the query not recognizing my column?


texts: Troubleshooting ArcObjects performance problems

texts:
We have developed a relatively complex geoprocessing routine using a C# / ArcObjects (COM) / ArcGIS Server 9.3.1 / ArcSDE 9.3.1 / Oracle 11g stack. Running the routine on a complete dataset takes many, many days and often does not complete.

Troubleshooting is hindered by the opaque nature of COM and the many levels of abstraction.

Specifics aside, could anyone suggest a heuristic for troubleshooting performance problems with ArcObjects? Hopefully the answer could be useful to others.

Edit: Using a profiler (JetBrains dotTrace), I can see that my problem calls are: ITable.GetRow, ITable.Select, ICursor.NextRow, ITable.Search, ISelectionSet.Search, ISelectionSet.get_IDs, IFeatureClass.Select, and IFeatureClass.GetFeature.


texts:
I have a 2017 Land Use raster of the US, shown below

that I'm trying to reclassify. It has 256 different values right now and my plan is to only have nine at the end - a few crops I care about and 'other'. I made a table in Excel that I'm trying to import into the reclassify tool so I won't have to change 256 values in ArcMap by hand.

The table is, for example,formatted like this:

And so on. When I try to load the table into the reclassify tool, it says "error, couldn't load table" and nothing else. I've tried it as a .xslx and a .txt file

Additionally, when I try loading the table to my ArcMap file first, it then becomes corrupted and I can't open that file from my desktop anymore.


texts: Cannot assign default value to layer created with processing tool - postProcessLayer

texts:
I am on QGIS 3.16 and I need to make a processing script that generate a new empty point shapefile with some fixed field, assigning them a default value and load it to canvas, so that if the user add a new features some fields are already compiled with the related values.
I have done several attempts but nothing seems to work... Is it possible?
This is my code
 def processAlgorithm(self, parameters, context, model_feedback):
    # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
    # overall progress through the model
    feedback = QgsProcessingMultiStepFeedback(1, model_feedback)
    results = {}
    outputs = {}

    layerFields = QgsFields()
    layerFields.append(QgsField('ID', QVariant.Int))
    layerFields.append(QgsField('Value', QVariant.Double))
    layerFields.append(QgsField('Name', QVariant.String))

    fn = os.path.join(parameters['Path'],parameters['Domain'] + '_sections.shp')
    writer = QgsVectorFileWriter(fn, 'UTF-8', layerFields, QgsWkbTypes.Point, QgsCoordinateReferenceSystem('EPSG:4326'), 'ESRI Shapefile')

    vl = QgsVectorLayer(fn, 'GUYGUY', 'ogr')
    idx = vl.fields().indexFromName('ID')
    vl.setDefaultValueDefinition(idx, QgsDefaultValue('$x', True))
    QgsProject.instance().addMapLayer(vl)

I slightly modified the approach and I understand I shomehow have to use the posprocessing capabilities... But I am not sure on how to deal with them... Are there any advice on how to start with them?
def processAlgorithm(self, parameters, context, model_feedback):
    # Use a multi-step feedback, so that individual child algorithm progress reports are adjusted for the
    # overall progress through the model
    feedback = QgsProcessingMultiStepFeedback(1, model_feedback)
    results = {}
    outputs = {}

    layerFields = QgsFields()
    layerFields.append(QgsField('ID', QVariant.Int))
    layerFields.append(QgsField('Value', QVariant.Double))
    layerFields.append(QgsField('Name', QVariant.String))

    fn = os.path.join(parameters['Path'],parameters['Domain'] + '_sections.shp')
    writer = QgsVectorFileWriter(fn, 'UTF-8', layerFields, QgsWkbTypes.Point, QgsCoordinateReferenceSystem('EPSG:4326'), 'ESRI Shapefile')

    context.addLayerToLoadOnCompletion(fn, QgsProcessingContext.LayerDetails('', QgsProject.instance(), ''))

    return results


texts:
can anyone please help with fitting together these 2 shape-files?
many many thanx

https://www.dropbox.com/s/saa9yp7ecg1qxh1/Grenzen.rar


texts: How to form SQL Expressions with mix of integer/numeric and text/string fields/values in ArcPy?

texts:
I am trying to create a query for the select by attribute function within a Python script.This is my current code.

    lot = '96'
    dp = 756421

    lotFile = "D:\\NCT_GIS\\DATA\\Regional_55.gdb\\Murray\\cadastre"
    arcpy.MakeFeatureLayer_management(lotFile, "lot_lyr")
    arcpy.SelectLayerByAttribute_management("lot_lyr", "NEW_SELECTION", "\"LOTNUMBER\" = " + lot + " AND " + "\"PLANNUMBER\" = " + str(dp))

Unfortunately I cannot concatenate and integer and if I convert DP to a string then the query is invalid as the PLANNUMBER field is numeric whilst LOTNUMBER is text.

Does anyone know how to get the variables into the query properly so that I can input any lot and dp as this will eventually become an addin where users can specify which land parcel they are interested in.


texts:
I have created a VRT using gdalbuildvrt from three single-band GeoTiffs (1-m resolution) as follows:

gdalbuildvrt chm_2157.vrt P_712806_CHM.tif P_712808_CHM.tif P_712810_CHM.tif  -a_srs 'epsg:2157' -srcnodata 1.175494351e-38

When I checked the output in QGIS, I noticed that a gap of 1 pixel (1-m) has appeared between two of the input tiles. Below are pictures of the input tiles (no gap present) and the VRT output (gap artifact included)

Input GeoTiff (no gap)

Output VRT (gap)

I've read threads about this issue occurring with gdalwarp and gdal_merge.py (e.g. here ) but the solutions don't appear to apply to gdalbuildvrt.


texts: I am having a small problem, hope everyone will help me! I want to create a map with a bar chart like in the demo of Amcharts 4 there is a map example with a pie chart. I don't know how to do so.this is the map I want to create, but a bar chart is not a pie chart


texts: Sort attribute table by hidden column in Print Composer

texts:
In a QGIS 2.18.20 attribute table I have some values with random numbers, but at the print composer I want to show them with order (clockwise). So I made a new integer column, and added numbers with the way of sorting I want to appear at the print composer.

Now I can't find a way to hide the column I used for sorting the numbers. I need to use it to sort the columns, but when I delete it or do something similar the order disappears. Is there a way to keep the sorting the way I want, but to hide the sorting column? Can I do it with a different way? In the picture, Field 1 is the point number, field2&3 are the coordinates and ORDER is the sorting column


texts: Calculating new field for raster table based on value from joined table using ArcPy?

texts:
I need to know whether it is possible to calculate a new field for a raster table based on a value from a joined table. Ideally I'd like to be able to accomplish this task without having to add table and raster layer to the table of contents within ArcMap. Please help!

The following section of my script runs successfully until it gets the line which creates the update cursor:

    #ArcMap v. 10.2.1
    # Add 3 fields to each of the new product rasters
    #Perform a spatial join connecting the variable rates table to the product raster

if sPotash and sPotash != "#":
    arcpy.PolygonToRaster_conversion(sPotash,"ZONE",""C:\Users\Desktop\CurrentWork\pm\Output_1192015\pt_tst","CELL_CENTER","ZONE",1)
    arcpy.AddField_management("pt_tst","LB_P_AC", "DOUBLE","","", "","")
    arcpy.AddField_management("pt_tst","LB_P_SQM", "DOUBLE","","", "","")
    arcpy.AddField_management("pt_tst","P_P_SQM", "DOUBLE","","", "","")

    arcpy.MakeRasterLayer_management(""C:\Users\Desktop\CurrentWork\pm\Output_1192015\pt_tst","pt_tst_lyr")

    arcpy.TableToTable_conversion("Variable_Rates.csv",""C:\Users\Desktop\CurrentWork\pm\Output_1142015","Potash_Rates","#","Potash")

    arcpy.AddJoin_management("pt_tst_lyr","Value","Potash_Rates.dbf","Zone")

    #arcpy.JoinField_management("pt_tst_lyr", "VALUE", "Potash_Rates.dbf","ZONE","Potash")

    potashrstr = pt_tst_lyr

    cursor = arcpy.UpdateCursor("pt_tst_lyr")

    fields = ["Potash","LB_P_SQM"] #the Potash field is a field belonging to the Variable_Rates table

    with arcpy.da.UpdateCursor("pt_tst_lyr",fields) as cursor:
        for row in cursor:
            row[1] = row[0]/4046.86


texts: How to find errors (gaps and overlaps) in a line feature class in ArcMap

texts:
I have a line feature class that contains roads that were manually edited. Now I want to check if all road crossings have been edited correctly. I want to check that there is no small gap or overlap (<1m) at the road crossings. Dead ends are not errors.
Which tool in arcmap could I use to highlight these errors?

To find small overlaps a solution is here, using first TOPOLOGY and second the tool "Feature Vertices to Points": Automatically marking all endpoints as exception in ArcGIS topology?

I am still looking for a soution for small gaps.


texts: Querying polygons that lie inside 4 points

texts:
I have 4 points that I always get, I would like to query if the polygon in the multipoint lies inside those 4 points. I’m using PostGIS, PostgreSQL.

I Have 4 Points which are dynamically updated I would like to query all the polgyons that lie inside those 4 points.

My Polgyon points are:

x = 16.07514381, double y = 47.95055008;
x = 16.24965286, y = 47.91411209;
x = 16.23201180, y = 47.82641983;
x = 16.14973259, y = 47.79749298;
x = 15.93255901, y = 47.89419937;

here is what I did so far

SELECT gid
FROM areas
WHERE st_contains(
         st_polygon(
            st_linefrommultipoint(
               areas.geom
            ),
            0
         ),
         st_mpointfromtext(
            'MULTIPOINT(1 1, 2 2, 3 3, 4 4)'
         )
      );


texts: Obtain monthly and yearly average values of EVI using MODIS data on Google Earth Engine

texts:
I am currently trying to export and obtain rasters of EVI for:

Every month of the year across 2000 to 2018 ie. If January has two MODIS EVI images, I would like to take the mean of the images for January and export one image for January 2000, and January 2001 and so on.
Similarly, I would like an averaged value of EVI for every month of the year across 2000 to 2018. In other words, a single EVI raster for January, a single EVI raster for February and so on.

I tried code provided in the answer here, but I am still unable to get a month by month as well as yearly average of EVI as needed above.

My code: https://code.earthengine.google.com/7c5c4eada7044b08d4c45342fd84d8d1


texts: Select / Deselect and Identify Features in QGIS

texts:
Sorry if this is all a bit 101, but I'm new to QGIS.

I basically have a large Arc project that I can access in QGIS for reference, so will only really be using it to identify features. For example to check the classification of a road line (which may be correct even if a lozenge on a map is incorrect, for example.)

I have my project open, but am unable to select any features from layers that are checked.

I only have the option to Deselect Features from All Layers. However, when I try to do that, nothing seems to happen. So I am unable to select any features I want to identify.

This must be something pretty simple, but some help would be much appreciated for a newbie!

Thank you.


texts:
I am making a code on C# to interpolate z values of a line on a DEM. Below is the code that I've written. I have also tried step by step debugging, the program works fine till the Geoprocessor execution. But on Geoprocessor execution, it crashes the ArcGIS. What could be the problem:::

IMxDocument CurrentDocument = (ESRI.ArcGIS.ArcMapUI.IMxDocument)ArcMap.Document;
            IMap MapSelected = CurrentDocument.FocusMap;

        //The code below is to fetch the Raster Layer having name equal to what selected in combo box
        int LayerIndex = 0;
        string LayerName;
        string comboText = (string)comboDEM.SelectedItem;

        for (int i = 0; i < MapSelected.LayerCount; i++)
        {
            LayerName = MapSelected.get_Layer(i).Name;

            if (LayerName == comboText)
            {
                LayerIndex = i;
                break;
            }
        }

        Geoprocessor GP = new Geoprocessor();
        ILayer RasterLayer = MapSelected.get_Layer(LayerIndex);

        InterpolateShape Interpolator = new InterpolateShape();
        Interpolator.in_surface = RasterLayer;

        //The code below is to fetch the Feature Layer having name equal to what is given in textbox
        string XS_Line = txtInputCrossSectionLine.Text;
        for (int i = 0; i < MapSelected.LayerCount; i++)
        {
            LayerName = MapSelected.get_Layer(i).Name;

            if (LayerName == XS_Line)
            {
                LayerIndex = i;
                break;
            }
        }

        ILayer FeatureLayer = MapSelected.get_Layer(LayerIndex);
        Interpolator.in_feature_class = FeatureLayer;
        Interpolator.out_feature_class = @"E:\For Arc GIS\Data\ResultFromTIN";

        GP.Execute(Interpolator, null);


texts:
How can I export the layout view elements of an ArcGIS 10.6 project for use in the QGIS 3.4.2 print layout? Apart from the shapefiles that accompany the project, I would like to be able to make the same map layout that I use in ArcGIS available for QGIS. At the moment it's just an .mxd file.

I've seen some plugins like SLYR (https://planet.qgis.org/planet/tag/slyr/) and Mxd2Qgs (https://github.com/oware/Mxd2Qgs) but don't know if these will accomplish what I need.

Would vector or image formats work, and can these layout paths help me?


texts: ST_Contains. Check for true/false in Python cursor

texts:
I would like to create a list of toponyms, which are within some distance from a point. I have calculated a buffer (dbrows_transf_back) and now I am checking if toponyms from a gazetteer (dbrows_os) are within this buffer. If I get "True" as a result, I will append the toponym to the name_list. The problem is that I don't know how to catch true and false. When I print dbrows_contain I get [(False,)] or [(True,)]. I have tried if dbrows_contains == 't': if dbrows_contains == '[(True,)]': and if dbrows_contains is True:
I am not sure, what else I can do. Any ideas?

name_list = []
for j in dbrows_os:
    name = j[0]
    toponym = j[1]
    sql_intersect = "SELECT ST_Contains(ST_SetSRID(ST_GeomFromEWKT('"+dbrows_transf_back+"'), 4326), ST_SetSRID(ST_GeomFromEWKT('"+toponym+"'),4326))"
    cur.execute(sql_intersect)
    dbrows_contains = cur.fetchall()
    conn.commit()
    if dbrows_contains is True:
        print "true"
        name_list.append(name)


texts: Incorrect geocode of state for sync tables

texts:
My sync table looks like this: http://i.imgur.com/ZfDp1L7.png

What should be Lebanon, Kanasas is being geocoded to the lat/lon of Lebanon, Kentucky.

This has happened with other cities that have the same name in multiple states.

How can I ensure that that the state column is properly read for the sync table?


texts: How to find nearest point for each point in current layer using ArcGIS Desktop?

texts:
I have a point layer and want to know the distance of nearest point of current layer for each point in the layer.

I found some tools in ArcToolbox about it , but they need two layers for input and they can find the distance of nearest object between that two layers. I tried to use them and i added my point layer in that tools ( Same layer for both inputs ) but they couldn't work properly because they found same point in second layer and gave me zero distance for every points.

How can i find distance of nearest point for each point in my layer ? It is better for me to find an answer in ArcToolbox . Something such as Tools or Models.


texts: Manage space in legend print composer QGIS

texts:
I am struggling to figure out how to manage space in the legend.

I have extra space between some items in my legend, which I would like to remove:

How can I manage these spaces?

I built my legend by adding an extra group layer with the legend items in it stored in 3 different layers (Infrastructure, BAL and Radians).

Every time I used a rule based classification, and it doesn't show any issues for the BAL and Radians groups.

But it doesn't work as I'd like in Infrastructure group. There is extra space between the two NRO and Infra layers.

How can I remove this extra space?

EDIT : here is a screenshot of the legend in the main QGIS window, displaying like I'd like in the print composer.


texts: How to restore the table window in ArcMap 10.2 (stuck in maximized view)

texts:
I've somehow managed to maximize the table window in ArcMap 10.2 (after choosing > Open Attribute Table from the TOC).

Whenever this window opens, it's maximized, and I can't see an option to restore it so it doesn't take up the full screen - the only option is to close it.

(A possible complication is that I'm running this in Windows 7 via Bootcamp on a Mac).


texts: How to calculate Fractal dimension from DEM

texts:
I have a digital elevation model of a sea floor topography, and I have already creat TIN, shp, ras of this zone (only X, Y, Z(depth))

To calculate the topography complexity, I think fractal dimension (2D and 3D)index is important to my research.

I use ArcGIS 10.1


texts: Migrating SVG-based styles for QGIS to other computers?

texts:
I'm working on a project that has involved creating many recolourable SVG icons and bringing them into QGIS. These are then used in layer symbologies, with these subsequent layer styles exported as .qmls and .slds.

When I open these .qmls on another computer the svg points are not attached - rather there is just a blank spot.

I've tried to get just the svgs to load by exporting the style library from the computer that all of the svgs were created on. However, when I imported this into another computers style library everything was once again blank. Putting the folder with all of the svgs in the C:\Program Files (x86)\Quantum GIS Lisboa\apps\qgis\svg directory is not working either...they show up in the option to select svg icon, however when a layers qml file is loaded it is once again blank.

How can I move these svg/style to other computers so that when I load a file with a set qml layer the proper svgs load as well?


texts: Close files handlers used by processing:runalg("qgis.difference"......) or os.popen("gdalwarp -q -cutline.......)

texts:
How can I do to close all files used by functions that I use in my plugin like

processing:runalg("qgis.difference"......)

os.popen("gdalwarp -q -cutline.......)

After runing my function, if I want to run it another time, I have to close QGIS and reload it again.


texts:
I am trying to create a code where a feature is created, and then automatically a popup window will open. Is there anyone who can assist me please?

Here is a part of the code and I want the popup to open right after that one.

var draw = new OpenLayers.Control.DrawFeature(
wfs_layer, OpenLayers.Handler.Point,
    {
        title: "Create a feature",
        displayClass: "olControlDrawFeaturePolygon",
        multi: true
    });

Later on I have a code for the popup:

//---POPUP fyrir Attribute---//
var select = new OpenLayers.Control.SelectFeature(wfs_layer,{hover: true, highlightOnly: true});
  map.addControl(select);
  select.activate();

    var selectClick = new OpenLayers.Control.SelectFeature(wfs_layer,{});
  map.addControl(selectClick);
  selectClick.activate();

  function onPopupClose(evt) {
  select.unselectAll();
}

wfs_layer.events.on({
featureselected: function(event) {
    var feature = event.feature;
        feature.popup = new OpenLayers.Popup.FramedCloud
            ("pop",
      feature.geometry.getBounds().getCenterLonLat(),
      null,
                '<form name="inputform">'+
                '

NAME:	\
                <input name="name" type="text" size="30" maxlength="30" value="'+feature.attributes.nafnfitju+'" >

'+
                '
Heimild:	\
                <input name="source" type="text" size="30" maxlength="30" value="'+feature.attributes.heimild+'" >

'+
                '
Code1:	\
                <input name="addrattara" type="text" size="30" maxlength="30" value="'+feature.attributes.addrattara+'" >

'+
                '
Code2:	\
                <input name="Source1" type="text" size="30" maxlength="30" value="'+feature.attributes.source1+'" >

'+
                '
Open/closed:	\
                <input name="open_closed" type="text" size="30" maxlength="30" value="'+feature.attributes.open_closed+'" >

</form>'+
                '
<input type="submit" value="Save" onClick="javascript:ourfunction(\''+feature.fid+'\')">
'+
                '
<input type="submit" value="Delete" onClick="javascript:delete(\''+feature.fid+'\')">
',
      null,
      true,
                onPopupClose
    );
    map.addPopup(feature.popup);
},

//destroy popup when feature is no longer selected. Prevents showing 2 Popups at the same time
featureunselected: function(event) {
var feature = event.feature;
        map.removePopup(feature.popup);
  feature.popup.destroy();
  feature.popup = null;
}
});

Hopefully someone can assist me, thanks in advance.

Sjonni, Iceland


texts: Using macports python 2.7 with macports qgis

texts:
I installed qgis through macports on my mac running OSX 10.8.2 but when I run the application (qgis 1.8) it says that I don't have python-gdal installed.

I searched google and it appears that it is because qgis is using the pre-installed version of python that comes with OSX. I have my own version of python 2.7 built through macports, and would like to use it with qgis if possible.

Does anybody knows how to tell qgis to use that version instead?

Thank you all for your consideration.


texts: Seeking list of countries updating their CRS like Australia has changed to GDA 2020

texts:
Australia has changed to GDA 2020.
Is there a site which lists whether any other countries are updating their CRS?


texts:
I used Google Maps as a basemap. In Map Composer window it is showing in the
correct location in QGIS, but when I try to export it, it is shifted.


texts:
I have 2 lines in 2 layer and I want to make a split where it is crossing each other/intersection point in Qgis 2.16

I want to make one line from the 2 line in the crossing point. I want the new line to be constant. I need the blue line in the pics after splitting. DONT know how to split it in the crossing point...

I have a good picture to see what I really want to do.


texts: How can I get my Qgis area calculations in metres rather than degrees? My shapefiles are projected using WGS 1984 UTM fuseau 35 Sud

texts: I am using Qgis 1.7.4 to calculate the area of polygons. My shapefiles have been projected using the coordinate system WGS 1984 UTM fuseau 35 Sud, but my area calculations are still displayed in degrees. how do I get them in metres?


texts: QGIS API Label expression for vector layer

texts:
I'm using QGIS 3.3 C++ API (python dev should be able to answer too as both C++ and Python API share the same elements names)

I'm trying to display a label on each point of the layer with this expression :

"'L: ' || \"Ref_track\" || ', PK: ' || \"PK\""

This expression works in QGIS 3.3 as

'L: ' || "Ref_track" || ', PK: ' || "PK"

I just adapted it for C++ string writting

here is how i'm doing it :

QgsPalLayerSettings mySettings;
QgsExpression myexp("'L: ' || \"Ref_track\" || ', PK: ' || \"PK\"");
mySettings.fieldName = myexp;

QgsTextFormat txtf;
txtf.setSizeUnit(QgsUnitTypes::RenderPixels);
txtf.setSize(20);
txtf.setColor("#ffffff");
mySettings.setFormat(txtf);

QgsVectorLayerSimpleLabeling* simpleL= new QgsVectorLayerSimpleLabeling(mySettings);

myPointsVectorLayer->setLabeling(simpleL);
myPointsVectorLayer->setLabelsEnabled(true);
myPointsVectorLayer->setCustomProperty("labeling/isExpression", true);

// commit changes
myPointsVectorLayer->commitChanges();

myPointsVectorLayer->updateExtents();
myPointsVectorLayer->triggerRepaint();

When i'm replacing my expression by "Ref_track" it works, but by "\"Ref_track\"" it doesn't ... WHY ?


texts: OGR: Load layer with sql statement and write to new data set

texts:
I am having difficulties after loading a shapefile with a sql statement with ogr.
The layer and its features seem to be correctly loaded. Still I am not able to write them into a new shapefile or pass the layer to the gdal_rasterize function.
The shapefile remains an empty file and the raster has the typical no input/something wrong happened values from 1.79..e+308 to -1.79..e+308 (not visualized in QGIS).
Code example for the creation of the shapefile (paths are omitted):

file = ogr.Open(file_path)
sql = 'Select * From {} where {} <= 5'.format(layer_name, "CS")
tmp_layer = file.ExecuteSQL(sql)
print(tmp_layer)
>>><osgeo.ogr.Layer; proxy of <Swig Object of type 'OGRLayerShadow *' at 0x00000214D99617E0> >
for feat in tmp_layer:
    print(feat.GetField('CS'))
    >>>(prints the field of the subset of features as expected)
drv = ogr.GetDriverByName("ESRI Shapefile")
out_ds = drv.CreateDataSource(output_folder + '/asd.shp')
proj = tmp_layer.GetSpatialRef()
layer = out_ds.CreateLayer('asd', proj, ogr.wkbMultiPolygon)
# copy the schema of the original shapefile to the destination shapefile
lyr_def = tmp_layer.GetLayerDefn()
for i in range(lyr_def.GetFieldCount()):
    layer.CreateField(lyr_def.GetFieldDefn(i))
#write features from tmp_layer to new layer
for feat in tmp_layer:
    layer.CreateFeature(feat)
out_ds = layer = None

What is my mistake?


texts: Python error opening plugin manager in Quantum GIS

texts:
I just attempted to update my version of QGIS from something like 2.2.0 (can't remember exactly what it used to be) to 2.6.0 and now encounter an error - copied below - when I try to open the Plugin Manager/Installer.

I uninstalled my previous version of QGIS first and, in addition to QGIS 2.6.0, downloaded the latest versions of Python, NumPy, matplotlib, and GDAL as recommended on the site.

I've tried reinstalling older versions of QGIS again, but the error pops up all the time now.

Sadly, I wanted to update the program in order to access the latest plugins.  I'm on Mac OS X 10.8.5.

An error has occured while executing Python code:

Traceback (most recent call last):
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/pyplugin_installer/installer_data.py", line 460, in xmlDownloaded
    if isCompatible(QGis.QGIS_VERSION, qgisMinimumVersion, qgisMaximumVersion):
  File "/Applications/QGIS.app/Contents/MacOS/../Resources/python/pyplugin_installer/version_compare.py", line 188, in isCompatible
    if len(maxVer)TypeError: object of type 'NoneType' has no len()

Python version:
2.7.9 (v2.7.9:648dcafa7e5f, Dec 10 2014, 10:10:46)
[GCC 4.2.1 (Apple Inc. build 5666) (dot 3)]

QGIS version:
2.6.0-Brighton Brighton,

Python path: ['/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins/processing', '/Applications/QGIS.app/Contents/MacOS/../Resources/python', u'/Users/Lucy/.qgis2/python', u'/Users/Lucy/.qgis2/python/plugins', '/Applications/QGIS.app/Contents/MacOS/../Resources/python/plugins', '/Library/Frameworks/SQLite3.framework/Versions/C/Python/2.7', '/Library/Frameworks/GEOS.framework/Versions/3/Python/2.7', '/Library/Python/2.7/site-packages/numpy-override', '/Library/Python/2.7/site-packages/matplotlib-override', '/Library/Frameworks/GDAL.framework/Versions/1.11/Python/2.7/site-packages', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python27.zip', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-darwin', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/plat-mac/lib-scriptpackages', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-tk', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-old', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/lib-dynload', '/System/Library/Frameworks/Python.framework/Versions/2.7/lib/python2.7/site-packages', '/Library/Python/2.7/site-packages', '/Users/Lucy/.qgis2/python/plugins/mmqgis/forms', '/Applications/QGIS.app/Contents/Resources/python/plugins/fTools/tools']


texts: Computing coordinates of point knowing coordinates of observer point, azimuth and distance to unknown point?

texts:
I have the longitude and latitude of a known point O(x,y) where x and y are obtained from transforming the longitude and latitude from EPSG 4326 to EPSG 32635. This is the observer point. I have another point A(x1, y1) whose coordinates I don't know, but I know the distance from O to A in meters and the azimuth from O to A in radians.

How can I compute the coordinates of A knowing only these?

A function in Qgis in C++ would be the most useful, or at least a mathematical formula would be great.


texts: Adding Lat and Long to attributed table in ArcGIS for Desktop?

texts: i am new ArcGis 10.2 user and I have the attributed data for the land cover type (I use Raster to point tool) but it does not have the exact latitude and Longitude at the point. Could anyone please help or suggest me some idea about it? I read the website, they said I have to use the VBA script to cope with this problem, I tried but it show .


texts: How to project and scale layer to cartesian (x,y) in QGIS?

texts:
I am new to QGIS and I have two layers I want to plot in QGIS but I somehow cannot convert one of them to other projection systems. This is a problem because each layer ends up in a different position and they should overlap:

Layer one: comes from a .csv and the polygons come in a WKT and are in cartesian format (long, lat). This is the layer that I cannot reproject to other systems.

EDIT I just checked and when loading the first layer, QGIS gives a message that says that this layer has no CRS. This confuses me even more...

Layer two: comes from a .shp file (with its correspondent .prj, .shx, etc.) and  this one I can convert to other projections but I can't get it to Layer One's projection coordinates

How do I solve this? The second layer I downloaded it from the spanish national statistics institute (www.ine.es/en/censos2011_datos/cen11_datos_resultados_seccen_en.htm) , the first one came as a given. Also if someone knows how to do it in R, answers are more than welcome.

Screenshots, note the coordinates that QGIS yields (first image is layer 1, second layer 2):

sample data:ID NOMBRE  PROVINCIA_INE   MUNICIPIO_INE   PROVINCIA_INE_ID    MUNICIPIO_INE_ID    SECCIONES_CENSALES  WKT
0-EU-ES-28-01-020-085   Miraflores de la Sierra Madrid  SOTO DEL REAL   28  28144   2802801001, 2803401001, 2804501016, 2806701001, 2806701002, 2806701003, 2808501001, 2808501002, 2808501003, 2812001001, 2814401001, 2814401002, 2814401003, 2814401004  POLYGON ((-3.799082 40.849834, -3.796926 40.848022, -3.792479 40.848271, -3.790517 40.846764, -3.787734 40.845321, -3.78545 40.838735, -3.774378 40.826942, -3.758486 40.820663, -3.744765 40.808853, -3.744925 40.807285, -3.740955 40.8054, -3.734056 40.805329, -3.732387 40.803447, -3.724855 40.801006, -3.723403 40.798152, -3.725922 40.793204, -3.73773 40.782672, -3.740098 40.777218, -3.737176 40.76905, -3.744816 40.76572, -3.736435 40.757916, -3.722519 40.740666, -3.715269 40.736091, -3.710177 40.729491, -3.717865 40.728543, -3.726083 40.73087, -3.7327 40.730485, -3.732848 40.735456, -3.738881 40.744419, -3.758577 40.758865, -3.759803 40.761586, -3.763947 40.761308, -3.771251 40.763863, -3.783859 40.770073, -3.802624 40.770618, -3.807198 40.777535, -3.811236 40.779251, -3.811895 40.785679, -3.813971 40.786912, -3.813479 40.796021, -3.828123 40.808769, -3.827749 40.814979, -3.830158 40.817935, -3.829713 40.830901, -3.815049 40.841148, -3.815182 40.84314, -3.808855 40.848724, -3.804239 40.849631, -3.802079 40.848544, -3.799082 40.849834))
0-EU-ES-28-01-011-086   El Molar    Madrid  EL MOLAR    28  28086   2800901001, 2805901001, 2808601001, 2808601002, 2808601003, 2810801001, 2810801002, 2812901004, 2814501002, 2816401001, 2816401002, 2816801001  POLYGON ((-3.593485 40.750705, -3.58733 40.748226, -3.582657 40.750591, -3.576659 40.750026, -3.572034 40.743977, -3.561801 40.743847, -3.55564 40.738086, -3.547488 40.734555, -3.54409 40.734615, -3.539326 40.73283, -3.530666 40.734047, -3.525766 40.732757, -3.528184 40.713665, -3.532074 40.711415, -3.544003 40.692216, -3.546959 40.691724, -3.553206 40.687127, -3.555312 40.682785, -3.553432 40.680278, -3.555185 40.668964, -3.554494 40.666867, -3.55023 40.664916, -3.556324 40.662384, -3.555827 40.660775, -3.555775 40.65033, -3.559066 40.650485, -3.566692 40.650369, -3.568148 40.643543, -3.577771 40.645808, -3.584164 40.645122, -3.584663 40.653534, -3.589333 40.655884, -3.593463 40.654804, -3.597052 40.65802, -3.601595 40.669762, -3.601294 40.673301, -3.604856 40.679733, -3.600745 40.696273, -3.60206 40.696628, -3.603361 40.710937, -3.607021 40.715053, -3.608081 40.71439, -3.605649 40.72052, -3.595097 40.725437, -3.591719 40.728969, -3.593485 40.750705))


texts: Change color of part of the line in QGIS

texts:
Imagine a few km long pipeline. Pipeline consists of 4 types of pipe. Those types are alternating.

0 - 200m  A type of pipe go. (Green color)
200 - 230m B type of pipe etc..(RED color)
230 - 380m C type of pipe etc..(YELLOW color)
.. 10 000m A type of pipe go. (Green color)

In QGIS I have whole line in one color, no point, just line*. Line is not
straight. Is it possible to change color, part by part?

*EDIT: its kml file so I have no data, just line


texts: Datum shift issue between OSGB1936 and WGS84

texts:
I have to work between OSGB1936 and WGS84 and BOTH of these data sets started out life as OSBG1936 BounaryLine data.

I thought I was just at the end of a full revision of my master mapping until I was presented with this

I use QGIS on both PC and Mac platforms and some of these files date back to work I did three years ago when I was much less proficient than the enthusiastic amateur I am now.

And whilst I have been through many many levels of processing and manipulation I rarely reproject.

For some reason back in the past I also seem to have invoked WGS 84 / UTM zone 30N ESPG:32630 (suspect that was me trying to do a best fit betwixt the two)

But I have OnTheFly projection turned on... this shouldn't happen surely?
And more importantly, How do I get them BACK to the same space?

I realise that I must have to do a form of re-projecting- but don't even know where to begin with determining what shits I need nor how to calculate them

Command line, GDAL, OGR and Python are all a scary closed book to me unless I have explicit instructions (I still don't know where or how I would enter a command line statement in QGIS!)- I'm a GUI man all the way.


texts:
I have created a MBTiles file from Tilemill and I am visualizing it through php MBTiles server. Its working fine. Now I want to count the loaded tiles and store it to a database. In a simple word I want to count how many tiles are loading when user is seeing the map. How can I do this?

Following is the php server code by which I am calling mbtiles

<?php
$zoom = $_GET['z'];
$column = $_GET['x'];
$row = $_GET['y'];
$db = $_GET['db'];
  try
  {
    // Open the database
    $conn = new PDO("sqlite:$db");

    // Query the tiles view and echo out the returned image
    $sql = "SELECT * FROM tiles WHERE zoom_level = $zoom AND tile_column = $column AND tile_row = $row";
    $q = $conn->prepare($sql);
    $q->execute();

    $q->bindColumn(1, $zoom_level);
    $q->bindColumn(2, $tile_column);
    $q->bindColumn(3, $tile_row);
    $q->bindColumn(4, $tile_data, PDO::PARAM_LOB);

    while($q->fetch())
    {
    header("Content-Type: image/png");
    echo $tile_data;
    }
  }
  catch(PDOException $e)
  {
    print 'Exception : '.$e->getMessage();
  }
?>

And using the following line to call tiles through that php file

var layer = L.tileLayer('../mbtiles.php?db=map.mbtiles&z={z}&x={x}&y={y}', {
                            tms: true,
                            maxZoom: 09,
                            minZoom: 05
                        });


texts: Interpolation problem: Spline with barriers

texts: I want to create a bathymetry map of a gulf using spline with barriers.When the interpolation finished the map in the gulf was created very well but a value is also created around the island. I tried to close the gulf with a line or making a polygon in order to "keep" the prices within it and failed again! Does anyone know the reason?


texts:
I have access to ArcGIS 10 through university, however my laptop runs Windows 8.  I have checked the requirements and they have listed the previous operating systems (excluding 8)

http://resources.arcgis.com/content/arcgisdesktop/10.0/arcgis-desktop-system-requirements

Will it run? If not, are there any ways round this (other than buying Arc 10.1 or installing an older operating system)?


texts: Change field names in SQL Server Management Studio

texts:
I am building a spatial view using SQL Server Management Studio 2014. One of the tables I am using to build this view has field names that are not compatible with esri; the field names have slashes, hyphens and spaces.

Is it possible to rename the fields of the view being created as the view is being created? So something like;

CREATE VIEW dbo.NewView AS
SELECT * FROM
dbo.tableA
INNER JOIN dbo.tableb
ON dbo.tableA.UID=dbo.tableB.UID
EXEC sp_rename 'dbo.NewView.Bad/field-name','GOOD_FIELD_NAME'

The EXEC sp_rename works if I run it by itself after the view has been created and on one field at a time, which isn't an ideal solution.

Or is there some way to run a find/replace on the field names of the view after it has been created?


texts: Dropping & inserting specific characters within a field being used for labeling

texts:
I would like to find out how to properly format a conditional expression in labeling, so as to only get specific characters back out of a field that has 9 characters in it, but don't wish to label them on the map in full.

Sadly, I am not the best in scripting languages, I know a bit SQL but I see labeling or does VBscript, Jscript, or Python.

I've thus far been able to drop the first two characters, but I am not familiar how to get to the part where I keep certain characters and skip others.

Example data---  [FACILITYID] WF2820001, WF3540034, etc...

What I have so far built in VBscript is as follows

Function FindLabel ( [FACILITYID] )
  FindLabel = mid ([FACILITYID],3,50)
End Function

This is at least dropping off the WF that doesn't need labeled. I wish to then have two characters remain, insert a '-' symbol, keep one character, drop the zeros and replace with a tilda ~ symbol and keep the last one or two numbers if their anything higher than 01.

Example data would have this input [FACILITYID] WF2820001, WF3540034 and would then come back with 28-2~1, 35-4~34.


texts: How to create a sequential unique id by groups?

texts:
I'd like to use the field calculator in QGIS to create a unique feature (polygon) id that reflects different levels in a geographic hierarchy.  So say I have a REGION attribute within which there are several AREAS (unique features).  I have alphanumeric codes for the REGIONS, but within each I want to number the AREAS in sequence from 1 to whatever.  The order doesn't matter.  So the final unique ids would look something like:

GY001,
GY002,
GY003...
DL001,
DL002...

If I didn't want to start over again for each REGION, I would just calculate it as

concat(REGION, $rownum+1)

But how to restart the numbering for each REGION?


texts:
I can open the layer on QGis but in MapInfo I got that error.

I already have other layers I can open without any problem, is just with one I get this error.


texts:
I am using ArcGIS Desktop 9.3 and WGS1984.

I have three municipalities with a total surface of 178 sq km. I want to analyze land cover of those municipalities using both LANDSAT images with a pixel size of 30 x 30 meters and SPOT with a pixel size 10 x 10 meters.

I have a raster of a whole county (> 1000 sq km) with land cover info derived both from LANDSAT and SPOT. When I clip the whole county raster with polygons of three municipalities, in case of clipping SPOT-based raster I get the correct surface of 178 sq km. That is the clipped raster area, the same as the surface of the three polygons/municipalities I am cutting with

If I clip the LANDSAT raster, the surface is 25% smaller, according to the attribute table! I know there should be differences because the pixel size is not the same and on the edges of the polygon I am cutting with some small discrepancies will appear.

But what happened with the 25% of the surface?

I really have no clue what to do, I have done such things before, but never had such different results.

Could it be something wrong with the data in the attribute table?


texts: ST_RELATE - Intersect ignore shared boundary

texts:
Using postgis I am trying to find polygons that intersect, ignoring polygons that only have a shared boundary (edge).

ST_Relate(g1, g2, 'T********') is a standard intersection.
ST_Relate(g1, g2, 'TF*F*****') standard intersection with boundary touching false

An image of what I am trying to achieve (from: http://daniel-azuma.com/articles/georails/part-6):

I am confident in the underlying geometries as they are created inside of postgis by taking the overlay similar to here.

Similar to: postgis st_relate polygon inside another and sharing a boundary though I want to ignore cases where the only topological connection is the boundary.


texts: How to order points on a road linestring to follow the road

texts: So I'm using a osm2pgsql schema. When I get a road, the points in the linestring don't have an order, so when I try to draw the road, it won't follow the road but will be a bunch of lines. Do I have to do the ordering myself? If so hhow? How do I order the points so they follow the road.


texts: How to depict line symbols with length from attribute

texts:
I want to depict angled drill holes as a horizontal traces (plan view).
I have a CSV file with X,Y,Azimuth and Horizontal length.

Using the simple marker I can get a line pointing in the right direction, starting from the X,Y coordinates. But even when I set the size to the Horizontal Length attribute in map units, the line are not the right length. I am using UTM21N/WGS 84, so I assume it is not a projection issue.


texts: Geodjango Exception when importing django.contrib.gis.gdal

texts:
I have installed python 3.7.5 64 bit. GDAL 2.4.4 and corresponding python gdal from http://www.gisinternals.com/query.html?content=filelist&file=release-1900-gdal-2-4-4-mapserver-7-4-3.zip
I have also set environment variables correctly. I did a gdalinfo --version from command prompt, i could get the version number. When i try to run django server using python manage.py runserver i am getting the below error

And in the command prompt i could see errors


texts:
I've been trying to export my geometries from GEE to drive (or any other method) by   converting my geometry to a FeatureCollection and using the following code

Export.FeatureCollection.toDrive({
  collection: geometry4sub,
  folder: "GEE data",
});

However, I get the following error message: Cannot read property 'toDrive' of undefined


texts: Allowing ArcPy script to continue if raster with no projection information encountered?

texts:
How do I allow an ArcPy script to continue if a raster with no projection information is encountered?

import arcpy, sys, os

InFolder = r"C:\test_in"
OutFolder = r"C:\test_out"
OutSR = arcpy.SpatialReference(26918) # NAD83 / UTM zone 18N

arcpy.env.workspace = InFolder

rasterList = arcpy.ListRasters()
for Ras in rasterList :
     inRaster = os.path.join(InFolder, Ras)
     spatial_ref = arcpy.Describe(inRaster).SpatialReference
     if spatial_ref.type <> "Unknown":
          outFileName = "{}.{}".format(os.path.splitext(Ras)[0], "jpg") # Remove .sid and add .tif to Ras filename
          outRaster = os.path.join(OutFolder, outFileName) # Path and Filename of output raster (.tif)
          arcpy.ProjectRaster_management (inRaster, outRaster, OutSR)

The code works now.


texts: Creating QGIS text label in composer calculating values from attribute table?

texts:
I am working with QGIS 2:14, for now, I have a qgs project with a lot of layers and I want to create a label in the print-composer with the result of a calculated expression from one of this layers.

In one of my layers (recycling_points.shp ) I have 150 records with the recycling points in the city. Then, for every record I have a field (num_containers) with the number of containers of every recycling point, and I want to create a text label in the print-composer, based on a calculated expression, with the number of containers in my city, but I don't know how can do it correctly.

I've tried with the next expression:

layer_property('recycling_points', count"num_containers")

but is invalid.

Any idea?


texts: Deleting vector layer from GeoPackage using PyQGIS

texts:
Programmatically I have added the vector layers to an existing GeoPackage in QGIS 3.2 using Python. Similarly, I would like to delete a layer from the existing GeoPackage using Python code. But there is no example/help found online.
Here is my code to add the selected layers to an existing Geo package:
layer = self.iface.activeLayer()
layer.startEditing()
gpkgPath='D:/AA/Layer/sample_gpkg'
options = QgsVectorFileWriter.SaveVectorOptions()
options.actionOnExistingFile =QgsVectorFileWriter.CreateOrOverwriteLayer
options.layerName = 'test'
_writer=QgsVectorFileWriter.writeAsVectorFormat(layer,gpkgPath,options)


texts: How to pass a parameter into Control.onAdd in leaflet

texts:
I have the following code:

var legend = new L.Control({position: 'bottomright'});
             legend.onAdd = function (map,max,min) {
             var div = L.DomUtil.create('div', 'info legend');
             var grades = [];
             console.log("Maximium ligne 139: "+ max);
             console.log("increament ligne 329" + min);
             var min;
             var ccount;
             var increament = (61544-4161)/8;
             for(var i=0;i<8;i++){
                 console.log("What we got! " +increament);
                 grades.push(Math.round(increament*i));
             }
            // grades = [0, 200, 1000, 1500, 2000, 3500, 5000],
            var labels = [];

                 // loop through our density intervals and generate a label with a colored square for each interval
                 for (var i = 0; i < grades.length; i++) {
                     div.innerHTML +=
                     ' ' +
                     grades[i] + (grades[i + 1] ? '&ndash;' + grades[i + 1] + '
' : '+');
                 }

                 return div;
             };

I'm trying to customize my map legend in leaflet. Therefore I need to pass some parameters and treat them inside the function in ligne 2.

I've tried to read the variables using this.max // this.min but its not working. Also I tried to pass the parameters inside the function (see line 2) but with no success. How can I achieve this ?


texts: How to load MrSID raster layer in QGIS on Mac OS?

texts:
I have just upgraded QGIS to 2.0 and now I can't load my SID raster layers.  I now remember that getting them into 1.8 was a performance but I can't remember all the details.

I have downloaded the latest SDK from Lizard -- I can remember that much and spent some time with google to try and figure out what to do next but I can't even find the 1.8 instructions.


texts: Uploading own areas 'videos' on Google Earth Street View?

texts:
In Google Earth you can drop a human icon to view a kind of video view of popular places. I guess these are basically converted from a video taken of the place.

If Google allows anyone besides itself to contribute content, like making a video of my own street and somehow make it available on Google Earth, then where can I read guidelines about how to do that?


texts: How to reclass aspect raster to float values

texts:
So I have an aspect raster that was created from a DEM. I'm trying to reclassify it using python, but I'm running into some issues. Here's the setup:

I need to reclassify values in the aspect raster that are between 315 and 360 to 1.2, between 0 and 45 to 1.2, between 135 and 225 to 0.8, and everything else to 1.0.

I've tried:

outAspect_Reclass = arcpy.sa.Reclassify(outAspect,"Value","0 44.99999 1.2;45 134.99999 1;135 224.99999 0.8;225 314.99999 1;315 360 1.2","NODATA")

And I've tried:

outAspect_Reclass = Float(Con((outAspect >= 0.0) and (outAspect <= 45.0), 1.2, Con((outAspect >= 135.0) and (outAspect <= 225.0), 0.8, Con((outAspect >= 315.0) and (outAspect <= 360.0), 1.2, 1.0))))

The first try fails because Reclassify won't reclassify to floats, and I'm not sure about the second try- can the Con() function do ranges (>=/<=)?

Can anyone tell me how to get this to work?


texts: GeoNetwork CSW Harvest Error - "gmd:MD_MetaData, csw: record" is not a valid value for QName

texts:
I'm trying to run a CSW Harvest in GeoNetwork, but I get this error:

"gmd:MD_MetaData, csw: record"  is not a valid value for QName

I'm using GeoNetwork 2.10.3 on Windows 7 64 bit and GeoServer 2.6.2


texts: search a point of interest using GeoServer

texts:
I develop an application of webmapping using GeoServer, GeoExt and ExtJS. I want to search my points of interests (layer that exist in GeoServer) so i create a two combox where i have to put the categorie of my points of interest and the name of my point of interest(optional). My data base is in PostgreSQL.

Any idea guys of how i should proced ?


texts: Problem PostGIS Spatial Query ST_Contains ST_Disjoint ST_Intersects

texts:
I'm working on a project with two different tables inside a database. One contains points in a column geometry(Point,4326) type and the second one has a column geometry(MultiPolygon,4326). Now I would like to do a query on pgAdmin III using st_contains and st_disjoint to understand if a point is inside or outside the entire polygon. Both tables have more than one row. The multipolygon table has inside values to represent the area, instead the point table has just points that are on the map.

Between tables there is no relation because the multipolygon table is created from a shapefile, instead the polygon table is created from collected data.

I tried to load both tables on QGIS and I can see that some points are over the area represented by multipolygon table.

I wrote this query:

SELECT distinct points.id, st_disjoint (points.geom, multipolygon.geom)
FROM points, multipolygon
ORDER BY points.id;

The result is this:

As you can see, the id 1,2,3,4 are ok and also checking on QGIS is fine. The number 5 and 6 are cointained inside the multypoligon, so the result would be false, but I get a false and a true result together.

Writing also:

SELECT distinct points.id
FROM points, multipolygon
WHERE st_disjoint (points.geom, multipolygon.geom)
ORDER BY points.id;

the result is wrong with all points inside the result. It's like the result I get with

WHERE st_intersect(points.geom, multipolygon.geom)='f'

If I write:

SELECT distinct points.id
FROM points, multipolygon
WHERE st_intersects(points.geom, multipolygon.geom)
ORDER BY points.id

The result is correct, without any repetition and it looks like I did something like the st_contains.

To do the st_disjoint, instead, I tried to do

WHERE not st_intersects(points.geom, multipolygon.geom)

or

WHERE st_intersects(points.geom, multipolygon.geom)='f'

but the result is the entire list of points. So it's not working to emulate st_disjoint.

Why the st_intersects sometimes works, when st_contains and st_disjoint don't work?


texts:
Here is my workflow for creating tiles and linking them to a layer.

1) load shapefile
2) grab layer representing shapefile
3) init & execute mapcruncher
4) turn on dynamic display
5) create dynamiccachelayermanager and init with my layer
6) save layer as layerfile

I've managed to successfully modify this sample:
http://help.arcgis.com/en/sdk/10.0/arcobjects_net/conceptualhelp/index.html#//000100000mm8000000

It creates my png files and a file called conf.xml.  In conf.xml it has the following data:

<LODInfo xsi:type="typens:LODInfo">
<LevelID>1</LevelID>
<Scale>485413.57628360949</Scale>
<Resolution>0.0011550226764674709</Resolution>
</LODInfo>

I thought if I set my mapScale to 485413.57628360949 then it would grab my generated tiles instead of regenerating them, but no luck, it still wants to regenerate them.  I think its just a mapscale issue, do you know how to get the mapscale of the generated tiles so I can put it set it via the IMapControl.mapScale interface?

Here's what happened next:

1) I used mapcruncher to generate tiles, it also generated a conf.xml
2) I open conf.xml to get the mapscales
3) I hard coded those mapscales in as my zoom levels
4) Dynamically generate the tiles at my hard coded zoom level
5) compare the mapcruncher generated tiles to the tiles which are dynamically created

The MapCruncher files start at Level 0 - 8, and the images are numbered in an orderly fashion.

C0000000.png C0000001.png... etc.

But the dynamically generated pngs are in the following format

C0000003d.png  even if only 1 png represents the whole map.   One would expect C00000000.png here....

I am truely lost here.


texts: I want to conect a slider with different columns of data in CartoDB

texts:
I am writing a map with CartoDB, and I would like to connect a slider with a select layer.
It is a slider of four points (2010-2011-2012-2013).
I want to connect the result of slider with the correspond column of database: data_2010; data_2011; data_2012 and data_2013
I am not sure if there is a function to do this, or if I have to use different layers in CartoDB, one to each data column.
This is the code:

    Time Slider | CartoDB.js

      html, body, #map {
        height: 100%;
        padding: 0;
        margin: 0;
      }
      #slider-range-max {
        position: absolute;
        bottom: 40px;
        right: 850px;
        left:110px;
      }

  #legend {
    font-family: serif;
    font-size: 27px;
    position: absolute;
    bottom: 100px;
    left: 100px;
  }
  #legend2 {
    font-family: serif;
    font-size: 27px;
    position: absolute;
    bottom: 70px;
    left: 100px;
  }

}

</style>
<link rel="stylesheet" href="http://libs.cartocdn.com/cartodb.js/v3/3.11/themes/css/cartodb.css" />
<link rel="stylesheet" href="http://code.jquery.com/ui/1.10.3/themes/smoothness/jquery-ui.css" />

    $(function() {
    $('#legend').html("Chile - Years")
    $('#legend2').html("" + 2010 + "     ---     " + 2011 + "     ---     " + 2012 + "     ---     " + 2013)
    $( "#slider-range-max" ).slider({
    range: "max",
    min: 2010,
    max: 2013,
    value: 2010,
    slide: function( event, ui ) {
    $( "#amount" ).val( ui.value );
    }
    });
    $( "#amount" ).val( $( "#slider-range-max" ).slider( "value" ) );
    });

<!-- include cartodb.js library -->
<script src="http://libs.cartocdn.com/cartodb.js/v3/3.11/cartodb.js"></script>
<script src="http://code.jquery.com/ui/1.10.3/jquery-ui.js"></script>

<script>
  function main() {
    var map = new L.Map('map', {
      center: [-39.66841891894786, -60.96533203125],
      zoom: 4
    });

    L.tileLayer('http://{s}.api.cartocdn.com/base-light/{z}/{x}/{y}.png', {
      attribution: 'CartoDB · Juan Moises de la Serna'
    }).addTo(map);

    cartodb.createLayer(map, 'http://juanmoises.cartodb.com/api/v2/viz/84a61b7e-b39c-11e4-bb07-0e853d047bba/viz.json', { legends: true })
     .addTo(map)
     .on('done', function(layer) {
    }).on('error', function() {
      cartodb.log.log("some error occurred");
    });
  }

  // you could use $(window).load(main);
  window.onload = main;
</script>


texts: GDAL's Gnomonic projection -- where do I set the elevation?

texts: I have an existing application that uses a Gnomonic projection where the elevation at the origin is added to the earth radius that is used in the projection calculations (6378137 plus elevation in m). I don't see anywhere to specify elevation/altitude in the SpatialReference.SetGnomonic method. Can I set some other attribute to achieve the same thing? Or can I modify the earth radius used by the projection?


texts: SSL certificate problem: unable to get local issuer certificate with gdal and curl

texts:
* Connected to ascendxyzdatastagingweu.blob.core.windows.net (168.61.57.78) port 443 (#0)
* libcurl is now using a weak random seed!
* SSL certificate problem: unable to get local issuer certificate
* Closing connection 0

I got the above eror with vsicurl and gdal and was not able to find any solution on google.

What is the tasks to enable gdal to use vsicurl with https?


texts: Union in QGIS: missing polygon boundaries

texts:
I am trying to UNION two shapefiles in qgis. I use the geoprossesing tool UNION. Both sf are in the same coordinate system (same as in the project).

In the resulting sf some polygon borders are missing (see mock-up of desired result in last image).

I also tried the overlay tool from Grass toolbox (or) but no result was produce!
Any ideas or a workaround ?


texts:
I try to export my map (layout or data view doesn't really matter) to a georeferenced tiff-file (geotiff). However, this presumably easy task is proving rather difficult to accomplish.

I first tried the code from this question:

import arcpy.mapping
mxd = arcpy.mapping.MapDocument("CURRENT")
arcpy.mapping.ListDataFrames(mxd)
arcpy.mapping.ExportToTIFF(mxd,r"D:\out2.tif","PAGE_LAYOUT",640,480,96)

All is good, a tif file is exported. But for some reason doesn't it contain any coordinates.

Some further research revealed this ArcGIS help page with a very detailed workflow. But it refers to an option that is not available in my export window:

  4, Choose Layout GeoTIFF (*.tif) from the Save as type drop-down list.

What am I missing here? Do I need another extension to create geotiffs?


texts: Change specific rows in attribute table using Arcpy in ArcGIS 10.3

texts:
I need to change values in field "name1" (string field) in specific rows: 70-72 into "no land use" value. I search Arcpy code to do it.

i red How to select specific row in attribute table using ArcPy (without iterating through whole table)? but didn't understand what shoud i do.

I using this code:

import arcpy

lyrName = r"G:\PROJECTS\herzel138\gis\layers\New File Geodatabase.gdb\digum2"
with arcpy.da.UpdateCursor(lyrName, ["name1"]) as cursor:
    for row in cursor:
        if row == 70 and row ==71 and row ==72:
            row[0] = 'no land use'
            cursor.updateRow(row)
            print "Processing complete"
        else:
            print 'nothing founded'

but i get no result:


texts: Creating points along 3D stream at inflection points

texts:
My goal is to generate points along my stream centerline at any positive inflection points (0.5 feet in elevation or greater). I used a Water Surface Grid (Raster) to convert my stream centerlines to 3D using the Interpolate Shape tool.
Looking at my stream in profile view, here are some examples of where I'd like to create points:

I posted a similar question before, but this script compares absolute difference between any two vertices and does not capture inflection points. Create points along 3D stream where elevation changes by 1ft
My idea is to store the starting value then look at the next vertex, if less than 0.5, skip and move onto the next vertex. Continue this till I hit a vertex that is 0.5 or greater from the starting point. Then create a point at this vertex. Then, move on and find the next jump of 0.5 or more.
Does someone have any ideas on how to script this out using ArcPy & Python 3?


texts:
I have a shape with 113 polygons. (partial catchment in hydrometric station sections). Is there any method, in ArcGIS, to overlap and merge polygons from upstream to downstream, in order to results a shape with 113 polygons not partial (between stations) but total (every partial catchment from downstream to be merged with the one upstream and so on....)

Below I selected a few partial catchments, as an example


texts: Multiple active TileLayers rather than just one?

texts:
I was wondering if it's possible to use multiple TileLayers in leaflet at the same time, rather than just having one active? A bit like group layers. In my map I have a base map, and I would like to make another TileLayer with transparent PNG:s where I can have roads, names of places, and so on.

Right now I am following this tutorial to achieve what I want:
https://www.hreikin.co.uk/2018/06/07/user-guide-how-to-create-an-interactive-campaign-map-with-leafletjs/

Basically I am working with the following tools:
Leaflet as base for the map
Zoomify (software) http://www.zoomify.com/ to creat the map Tiles from the image I have made in Photoshop.

The base tile looks like this in the code:

var layer = L.tileLayer.zoomify('./basemap/{g}/{z}-{x}-{y}.jpg', {
    width: 19000,    // MUST be defined.
    height: 10688,   // MUST be defined.
    minZoom: 2,
    maxZoom: 7,
    continuousworld: false,
    noWrap: true,
}).addTo(map);

Then I have a bunch of markers for cities that are grouped into layers. What I currently would want to do is to make multiple basemaps as layers, rather than toggleable.

As you can see on the layer icon in the top right corner you can only have one active basemap at a time: https://hreikin.co.uk/interactive-campaign-map/interactive-campaign-map.html

Does anyone have any suggestions?


texts: Calculate how to cover an area by recepetion tovwers

texts:
I am working in a Nature Reserve and we are putting up towers that will be used to to track certain animal with collars on.
The towers work on line of site only. We must place the towers in a matter with the least overlap and cover the whole park area. More smaller areas will have the most accurate and current results.

Things to consider are:

Terrain: will break line of site
Hight of tower - will determine area covered

How would you go about creating reception based on line of site and the terrain?

I understand first we would need some kind of terrain map with contours etc. The problem is the ones I can find are only 5m intervals, but a 1m change in tower hight would have a huge impact on the coverage of the tower.

We must assume power settings and range are theoretical, and will be tested in the field


texts: pyqgis: displaying 2 different datasets in one instance

texts: I have a need to display 2 completely different sets of spatial data in a single standalone qgis application. This means, 2 separate sets of layers, 2 map canvases, etc. Is this at all possible? (since some of the core objects in a qgis application are singletons, ie. QgsMapLayerRegistry)


texts: AttributeError from PDF Export from ArcMap Bookmark Loop

texts:
import arcpy, os, datetime

folderPath = r"C:\path"
dateTimeStamp = time.strftime('%Y%m%d')

basename = "Flows"
mxd = arcpy.mapping.MapDocument(r"C:\path\file.mxd")
df = arcpy.mapping.ListDataFrames(mxd, "Layers")[0]
#loops through bookmarked locations and executes a PDF export
for bkmk in arcpy.mapping.ListBookmarks(mxd, "",data_frame=df):
    #going through text elements to find the map title and set it to the bookmark name
    for elm in arcpy.mapping.ListLayoutElements(mxd, "TEXT_ELEMENT", "documentTitle"):
        elm.text = str(bkmk.name) + " - Flows"
        #set extent to bookmark extent before pdf export
        df.extent = bkmk.extent

    print basename + "_" + bkmk.name + " is exporting to PDF"
    arcpy.mapping.ExportToPDF(mxd, basename + "_" + bkmk.name + "_" + dateTimeStamp)
    print basename + "_" + bkmk.name + " successfully exported to PDF"
    print "Done exporting this folder to PDF - " +folderPath

#cleanup
del mxd

I'm receiving the following error after attempting to loop through about 10 bookmarks worth of PDF exports. The odd thing is, if I do 3-4 bookmarks at a time, it will export successfully but not if I do the full set of 10.

  AttributeError: PageLayoutObject: Error in executing ExportToPDF


texts: How to get countries through OpenLayers.Protocol.WFS?

texts:
Am able to get states information of United States by using OpenLayers Protocol.

new OpenLayers.Protocol.WFS({
                   version: "1.1.0",
                   srsName: "EPSG:900913",
                   url:  "http://demo.opengeo.org/geoserver/wfs",
                   featureType: "states",
                   featureNS: "http://www.openplans.org/topp"
               })

I need to get Countries information, tried replacing the fetureType: "countries", but this dint work. Can anybody help me.


texts: Where to install python plugins without internet access?

texts: I have a problem where I have no internet access. In the thread regarding a similar subject, it was mentioned that the plugin .zip file should be extracted to the ".qgis/python/plugins" in the home directory. Where would I find the "home" directory.


texts:
i know in this site everybody more or less have been setup a geoportal for showing their basemaps and analysis on local or net for their workers or conserned people in the net.

it is obvious that in the gis platform there are lots of datas and machines for serving all requests from people. if you are alone in your work or more knowledgeable from other workers in the same place it is meant to there are lots of works to do in your professional businness if you respect your work...

owing to you loving your work, doing more work sth like coding, designing or mapping, working a lot probably not bother you. but when you are sleeping at middle of the night or with your favorite friends come together with stealthy way for doing sth. , so at your special moment when comes a call about job makes you bother or it disgusts you from your work...

at now for managing and controling the gis sth like only restarting services as apache, nginx or iis; your database maybe swelled; maybe image server not working bla bla bla... what tools are you advice me for controling not everything but more in gis or for gis....do you have any open source tool for this..

a program with defining database controlling size of data, service status or restarting them...

thank you for your advice....
thanks.


texts: How to find Coordinate System of datasets inside MDB and GDB

texts: I have been given some ESRI personal geodatabase files (MDB) and file geodatabase files (GDB). I don't know the coordinate system(s) for these files. Is there a way to do this without using ESRI software?


texts: How save modified features in server side with java and geotools?

texts:
I'm developing a gis Application with java. I use Spring 4.0.1 and GeoTools-13-RC1 in server  side and OpenLayers 2.13.1 in client side. In client side, I have a wfs layer and it is as follow:

var saveStrategy = new OpenLayers.Strategy.Save();

wfs_Layer = new OpenLayers.Layer.Vector("Editable Features", {
    strategies: [new OpenLayers.Strategy.BBOX(), saveStrategy],
    projection: new OpenLayers.Projection("EPSG:4326"),
    protocol: new OpenLayers.Protocol.WFS({
        version: "1.1.0",
        srsName: "EPSG:4326",
        url: 'point/read',
        featureType: "point",
        featureNS: "http://itsme.opengeo.org",
        geometryName: "the_geom",
        readFormat: new OpenLayers.Format.GeoJSON({}),
    })
});

wfs_layer get all features from server. Sever return features with follow code:

try {
    Map<String, Object> params = new HashMap<String, Object>();
    params.put("dbtype", "postgis");
    params.put("host", "localhost");
    params.put("port", 5432);
    params.put("schema", "public");
    params.put("database", "template_postgis");
    params.put("user", "postgres");
    params.put("passwd", "postgres");

    DataStore dataStore = DataStoreFinder.getDataStore(params);

    SimpleFeatureStore pointStore = (SimpleFeatureStore) dataStore.getFeatureSource("point");

    SimpleFeatureCollection featuresCollection = pointStore.getFeatures();
    SimpleFeatureIterator features = featuresCollection.features();

    FeatureJSON fjson = new FeatureJSON();
    StringWriter writer = new StringWriter();
    fjson.writeFeatureCollection(featuresCollection, writer);

    String json = writer.toString();
    return json;

} catch (IOException ex) {
    Logger.getLogger(PointController.class.getName()).log(Level.SEVERE, null, ex);
}

The returned features format is as bellow:

{
    "type": "FeatureCollection",
    "features": [{
            "type": "Feature",
            "geometry": {
                "type": "Point",
                "coordinates": [40, 45]
            },
            "properties": {
                "name": "itsme 1"
            },
            "id": "point.1"
        }, {
            "type": "Feature",
            "geometry": {
                "type": "Point",
                "coordinates": [35, 40]
            },
            "properties": {
                "name": "itsme 2"
            },
            "id": "point.2"
        }]
}

User can edit this layer in client side and with saveStrategy.save() send edited features to server. Browser send edited features as follows:
Request Payload:

<wfs:Transaction xmlns:wfs="http://www.opengis.net/wfs" service="WFS" version="1.1.0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.opengis.net/wfs http://schemas.opengis.net/wfs/1.1.0/wfs.xsd">
    <wfs:Update typeName="feature:point" xmlns:feature="http://iran.opengeo.org">
        <wfs:Property>
            <wfs:Name>the_geom</wfs:Name>
            <wfs:Value>
                <gml:Point xmlns:gml="http://www.opengis.net/gml" srsName="EPSG:4326">
                    <gml:pos>48.694467049643556 35.41594833365911</gml:pos>
                </gml:Point>
            </wfs:Value>
        </wfs:Property>
        <wfs:Property>
            <wfs:Name>name</wfs:Name>
            <wfs:Value>itsme edited</wfs:Value>
        </wfs:Property>
        <ogc:Filter xmlns:ogc="http://www.opengis.net/ogc">
            <ogc:FeatureId fid="point.2"/>
        </ogc:Filter>
    </wfs:Update>
</wfs:Transaction>

Now I want to modify featuers in server side. How do I do?


texts:
Im using leaflet to draw points defined by custom div(like the default css in the leaflet.css).
I have to update the points every 30 mili secs.
But it seems to be so slow.. when i reach to high number of points.

I thought about forcing leaflet to use the gpu somehow instead of using the cpu to render the points.
Try to use transalte etc.. but it doesnt show any changes.

So there is anyway of using the gpu to render those points and makes my app to work faster?


texts: Exporting layer symbolized in two different ways with qgis2web?

texts:
I have a layer containing all the counties of the US. I would like to symbolize this layer for presentation based on two different attributes, say total population and population per square mile, using colour ramps for each.

Only one symbolization would be presented at a time, so one would view one symbolization or the other by turning on or off the higher layer.

I want to export this map using qgis2web. I can achieve this effect by duplicating the layer and setting the symbolization appropriately, but this results in two large .js-files being written containing the feature sets, which are identical.

Is there a way to symbolize a single layer to achieve this effect? Rule-based symbolization does not achieve this goal as far as I can tell.

I have QGIS 3.4.2 and the latest version of qgis2web


texts:
I'm not very used to JavaScript, and I don't understand how to adapt my code to this tutorial : https://docs.mapbox.com/mapbox-gl-js/example/popup-on-hover/

Indeed, in my code, the text information is from a layer made in Mapbox studio :


texts: QGIS and Google map doesn't match/lineup!

texts:
I am trying to draw a map but in QGIS the map doesn't correctly show one of my study area (Station-4) correctly; however, the Google Map does.

What should I do now to get a correct map (Station-4) just like the one shown by Google Map.


texts: How to create regular point grid inside a polygon in Postgis without knowing cell size but cell number?

texts:
I need to create a grid within a polygon like it is done in

How to create regular point grid inside a polygon in Postgis?

There a function is created to draw the grid points. But necessary to know is how huge the cells are or how huge is the distance between the point. In my case I don't know anything about that. I just know how much grid cells should be created. How can I solve that problem? How is it possible to create a grid without knowing the cell size but the number of cells that should be created? Does anybody can help me?


texts:
I've written a tool to create shapefiles using GDAL and am viewing the results in QGIS. Each time the tool runs, I want to create a new output directory*. I'm using QGIS 2.18.9 on Windows 10.

The problem is that, if I have viewed the previous result in QGIS, I cannot delete the output directory unless I close the whole application. Just removing the layer from the display isn't enough - there is still a lock on the directory. I don't want to close QGIS each time I run the tool. The browser panel doesn't seem to have a delete option, like ArcCatalog does, unless I'm missing something.

I have found that if I right click on the output directory and select Hide From Browser, I can run the tool again, but I don't know how to unhide it again. I have to click Add Layer and browse to the shapefile if I want to look at the output.

So the question is, how do I delete a folder, having viewed the contents in QGIS? I've used ArcGIS a lot, but I'm relatively new to QGIS. Is there a way to either delete the folder from within the browser or otherwise release it so I can delete it from elsewhere?

*This is by design - I want to put the onus on the user to make sure the output directory doesn't exist, rather than have the tool potentially delete existing data.


texts:
Is it possible to generate an API-key in Cartodb with read-only permission?

Seems like the default API-key has full controll


texts: Change floors on google maps street layer

texts: I'm using Web > OpenLayers plugin > Google Maps > Google Street to pull in the floorplan a user submitted to Google on a building.  I'm using it to trace my own shapefile.  But, it only shows the first floor.  On Google Maps web application, you can change the floor when you have the building selected, is there an option to do this in QGIS?


texts: PYQGIS signal for layer tree when active layer changed

texts:
I want to know if there is signals in layers tree to get notified when the active layer changed or more precisely when the user clicks or dbl click on any layer or node in the layers tree view.

I search in QGIS documentation but found nothing!


texts: GPS and directional layer for GPS points

texts:
I have millions of GPS location data for some vehicles and a street network for a city.

PROBLEM:
In an intersection approach, there are points which are entering the intersections and there are other set of points which are leaving the intersection. I want to draw a polygon and clip only those points which are entering the intersection. How this can be done ?
I am thinking of a directional layer which can only capture the points moving in one direction  ?


texts: Reproject all layers in a spatialite file using ogr2ogr

texts:
I would like to reproject all layers in a spatialite file from one srs to another srs.
Is there a way I can archieve this without parsing the layers names to a txt-file?
I know I can write all layernames to a txt file like this:

ogrinfo myfile.sqlite >layerlist.txt

But the content of the textfile has to be edited manually before I can use it to loop over the layernames as there is unneccessary content like this:

INFO: Open of `myfile.sqlite'
      using driver `SQLite' successful.
1: layername1 any (Multi Line String)
2: layername2 any (Multi Polygon)

Is there a way I can loop over the layers without having to edit the layerlist manually?


texts: Planet not able to get a response from quick-search API call

texts:
My question is related to:
Planet API not being accepted for quick search but is on stats and planet quick-search api not recognizing authentication?.

I am not able to successfully get a response from Planet 14 days trial API call for quick-search, I get the following error:

  {'message': 'Please enter your API key, or email and password.', 'errors': []}

Edit1: my code is as follows:

 #AREA of Interest
geojson_geometry = {
  "type": "Polygon",
  "coordinates": [
    [
      [-19.423828125,-4.653079918274038],
      [51.591796875,-4.653079918274038],
      [51.591796875,45.336701909968134],
      [-19.423828125,45.336701909968134],
      [-19.423828125,-4.653079918274038]
    ]
  ]
}
#Images overlapping out area of interest
> geometry_filter = {
>     "type": "GeometryFilter",
>     "field_name": "geometry",
>     "config": geojson_geometry }

#for a specific date range

> date_range_filter = {
>     "type": "DateRangeFilter",
>     "field_name": "acquired",
>     "config": {
>         "gte": "2009-01-01T00:00:00.000Z",
>         "lte": "2019-03-31T00:00:00.000Z"
>     } }

#images with less than 50% cloud range

> cloud_cover_filter = {
>     "type": "RangeFilter",
>     "field_name": "cloud_cover",
>     "config": {
>         "lte": 0.5
>     } }

#combine geometry, date and cloud filters

> combined_filter = {
>     "type": "AndFilter",
>     "config": [geometry_filter, date_range_filter, cloud_cover_filter] }
>
> import os import json import requests from requests.auth import
> HTTPBasicAuth

#PLANET_API_KEY = os.getenv('c43b6ae27338421dad81dd7c2a5306a6')

> os.environ["PLANET_API_KEY"] = "c43b6ae27338421dad81dd7c2a5306a6"
>
> item_type = "PSScene4Band"

#API request object

> search_request = {
>     "interval": "day",
>     "item_types": [item_type],
>     "filter": combined_filter } search_result = requests.post('https://api.planet.com/data/v1/quick-search',
> auth=HTTPBasicAuth(PLANET_API_KEY, ''), json= search_request)
>
> print(json.dumps(search_result.json(), indent=1))


texts: How to position contour labels in Maplex above the line?

texts:
I would like to utilize the "contour placement" options of "Uphill alignment" and "Laddering" while positioning the label above the contour not over it. This does not seem to be an option under the contour placement where it's forcing me to place the contour over the line. (Which makes no sense in ArcGIS since there is no easy way to mask that portion of the line efficiently. Possible but very cumbersome and not dynamic)
I tried the other placement options in Maplex but then "Uphill Alingment" and "Laddering" is not longer an option.

Does anyone has a good way of doing this?


texts: converting polygon data points (shapefile) into point datapoints (another shapefile), using QGIS

texts:
I am creating my own dataset for landuse classification.
I have four different shapefiles for 4 different categories, forest, water, park and buildings.
All these shapefiles contain polygon geometries of various sizes of respective categories.

I want to create shapefile of each category which contains points instead of polygons.
For example, my forest shapefile contains 10 polygons and i want to create another shapefile for the forest, which will contain 100 or more points. Obviously these points should lie within these polygons.

I don't want to manually select the points which will be quite cumbersome given the no of shapes and points.

I would like to know any  automatic method using QGIS.


texts: Extracting values from raster according to Lat and long of values?

texts:
having 12 files with 12 hdrs for one year:    these files are raster (projected  WGS84,lat long):

          samples = 1440
           lines   = 720
           bands   = 1
           header offset = 0
           file type = ENVI Standard
           data type = 4
           interleave = bsq
            byte order = 0
      map info = {  Geographic Lat/Lon, 1, 1, -180, 90, 0.25, 0.25,WGS-84}
        coordinate system string = GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",
                 SPHEROID["WGS_1984",6378137,298.257223563]]
       ,PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]]
                 }

These lines will open the files as a list:

      a<-list.files("D:\\ECV\\2010", "*.envi", full.names = TRUE)
           for(i in 1:length(a)){
            d <- raster(a[i]}

I would like to extract the values correspond to  44.8386° N, 0.5783° W from all files as txt file


texts: Checking if optional parameter/field has been used or not in ArcPy?

texts:
I have 6 arcpy.GetParameterAsText

The first one is a feature class and the remaining 5 are fields (some being optional).

I plan to use an if statement to check if any fields have been selected for the optional parameters.

I'm new to if statements in arcpy and have had a google. To test if things work I have:

if FWHier == '':
    print("Its blank")

There are no complaints about the syntax here but it does not return (print) anything.

I'm also wondering how to 'end' and if statement. Equivalent to the END IF in Mapbasic.


texts: QGIS equivalent to ArcGIS Aggregate Polygons?

texts:
I am looking for a tool in QGIS 2.40 that performs the function ArcGIS refers to as Aggregate Polygons. Basically it takes a group of polygons within a specified distance and creates a concave hull around them.

Does anyone know if this functionality (or anything similar) exists within QGIS currently?


texts: CQL_FILTER on multilayer request where only one layer needs filtered

texts:
I have an image mosaic timeseries layer called agdd as well as a shapefile layer called cb_2014_us_state_500k that displays us state boundaries.

I can stack both layers with the following query:

http://server/geoserver/gdd/wms?service=WMS&version=1.1.0&request=GetMap&layers=gdd:agdd,gdd:cb_2014_us_state_500k&time=2016-01-10&styles=&bbox=-125.020833333333,24.0625,-66.479166666662,49.937500000002&width=1400&height=700&srs=EPSG:4269&format=application/openlayers

The result looks like so:

I can use CQL_FITLER=NAME='Texas' to only show the Texas border in the cb_2014_us_state_500k layer:

http://server/geoserver/gdd/wms?service=WMS&version=1.1.0&request=GetMap&layers=gdd:cb_2014_us_state_500k&CQL_FILTER=NAME=%27Texas%27&styles=&bbox=-179.148909,-14.548699,179.77847,71.365162&width=768&height=330&srs=EPSG:4269&format=application/openlayers

Which looks like this:

Now I'm trying to both the combine layers AND filter the states layer to Texas all in one request. I tried to use a semicolon in the CQL_FILTER to specify that there are two layers, but it doesn't parse:

http://server/geoserver/gdd/wms?service=WMS&version=1.1.0&request=GetMap&layers=gdd:agdd,gdd:cb_2014_us_state_500k&CQL_FILTER=;NAME='Texas'&time=2016-01-10&styles=&bbox=-125.020833333333,24.0625,-66.479166666662,49.937500000002&width=1400&height=700&srs=EPSG:4269&format=application/openlayers

Any ideas how to write this request?

Here is the error I get when using an empty filter for one of the layers:

<?xml version="1.0" encoding="UTF-8" standalone="no"?><!DOCTYPE ServiceExceptionReport SYSTEM "http://geoserver-dev.usanpn.org:80/geoserver/schemas/wms/1.1.1/WMS_exception_1_1_1.dtd"> <ServiceExceptionReport version="1.1.1" >   <ServiceException>
      Could not parse CQL filter list.
Encountered &quot;&lt;EOF&gt;&quot; at line 1, column 12.
Was expecting one of:
    &quot;not&quot; ...
    &quot;id&quot; ...
    &quot;in&quot; ...
    &lt;IDENTIFIER&gt; ...
    &lt;DATE_TIME&gt; ...
    &quot;(&quot; ...
    &quot;[&quot; ...
    &lt;DATE&gt; ...
    &quot;-&quot; ...
    &lt;INTEGER_LITERAL&gt; ...
    &lt;FLOATING_LITERAL&gt; ...
    &lt;STRING_LITERAL&gt; ...
    &quot;true&quot; ...
    &quot;false&quot; ...
    &quot;point&quot; ...
    &quot;linestring&quot; ...
    &quot;polygon&quot; ...
    &quot;multipoint&quot; ...
    &quot;multilinestring&quot; ...
    &quot;multipolygon&quot; ...
    &quot;geometrycollection&quot; ...
    &quot;envelope&quot; ...
    &quot;include&quot; ...
    &quot;exclude&quot; ...
     Parsing : NAME=&apos;Ohio&apos;;.
</ServiceException></ServiceExceptionReport>


texts: gdal/ogr python adding new feature into existed PostGIS table

texts:
I am trying to add data into Postgis via gdal/ogr in python. However, if I first use OGRDataSource.CreateLayer("test") to create layer, then I can add new feature sucessfully. However, when tables are already defined in the database, if I use OGRDataSource.GetLayer(table_name) to get the layer and then try to add/edit feature. It does not work.

I do found differences after running OGRLayer.TestCapabilities() on both layers. If layer is reference by CreateLayer(),  the capabilities such as SequentialWrite and RandomWrite are both True. If the layer is reference to existed tables using GetLayer(), at least those two capabilities are False. Which I think explained the reason why add feature fails. But how can I enable them or what is the better way to insert new record(no geometry) into Postgis in Python?

**# The user is the owner of this database**.
conn_str = "PG: host=%s dbname=%s user=%s password=%s" %('localhost', 'CEdB', 'test', 'test')
conn = ogr.Open(conn_str)

def this_create_layer_and_add_feature_OK():
    # here we create test_made table
    layer = conn.CreateLayer('test_made', None, ogr.wkbUnknown, ['OVERWRITE=YES'])
    # add new fields
    err_code = layer.CreateField(ogr.FieldDefn('object_id', ogr.OFTInteger))
    field_def = ogr.FieldDefn('name', ogr.OFTString)
    field_def.SetWidth(50)
    err_code = layer.CreateField(field_def)
    feature = ogr.Feature(layer.GetLayerDefn())
    feature.SetField('object_id', 1)
    feature.SetField('name', 'testing')
    layer.StartTransaction()
    err_code = layer.CreateFeature(feature)
    feature = None
    layer.CommitTransaction()

def this_not_create_feature_in_table():
    # test_table is a predefined table, with 2 columns, no geometry!
    lyr = conn.GetLayer("test")
    new_feature = ogr.Feature(lyr.GetLayerDefn())
    new_feature.SetField("object_id", 1)
    new_feature.SetField("name", 'testing')
    lyr.StartTransaction()
    if lyr.CreateFeature(new_feature) == ogr.OGRERR_NONE:
        new_feature = None
        err_code = lyr.CommitTransaction()


texts: Scaling issues when exporting FBX models from CityEngine

texts:
I have extruded a dataset of circles into cylinders based on values in an attribute table using CityEngine and I am trying to export this into an FBX model for use in 3DS Max as well as Infraworks 360. The model is correct in CityEngine, however when I export the model the units and therefore scale are well off. The units of the model in CityEngine are in metres as it is set in the British National Grid Co-ordinate system, but when the FBX model is brought into 3DS Max it is being read as centimetres and converting it to metres then throws off the co-ordinates by factor of 100 so it is not properly geo-located.

I have read some forums online which suggest that CityEngine can only export in centimetres, does anyone know if this is true?

I have attached screenshots of the model in CityEngine, export settings for the FBX model and the scale issues in 3DS Max.


texts:
I am writing an addin, which needs to split a polyline into a number of segments taking the number of segments as an input from user.
How can I accomplish this?


texts: cartodb static named map with parameters

texts:
I have a named map with a variable that allow me to change the css_color value of a layer when I render it using cartodb.js.

I would like to create a static map for this named map but I can't find how pass the color parameter using the map API.

https://myusername.cartodb.com/api/v1/map/static/named/my_named_map/800/600.jpg
gives me :

{
  "error": "Invalid css_color value for template parameter 'color_up': "
}

https://myusername.cartodb.com/api/v1/map/static/named/my_named_map/800/600.jpg?config={color_up:"#f2bb00"}
gives me :

{
  "error": "malformed config parameter, should be a valid JSON"
}

I also tried ?color_up="#f2bb00" without success. Any ideas ? There is nothing in the documentation so it may not be possible to pass a parameter at all ?


texts:
I have a UN-Adjusted Population Count dataset, for the year 2015, stored in a TIFF file format (~450MB) downloaded from http://sedac.ciesin.columbia.edu/data/collection/gpw-v4/sets/browse

In R I am able to read the data using raster package. I am also able to plot it as follows,

#read your TIFF file

worldpopcount <- raster("gpw-v4-population-count-adjusted-to-2015-unwpp-country-totals_2015.tif")

worldpopcount

plot(worldpopcount)

Now I want to crop the West African region (covering the three countries Guinea, Liberia and Sierra Leone). My coordinates for West Africa are,

# extent format (xmin,xmax,ymin,ymax)

# WestBound     -16
# EastBound      -7.25
# SouthBound      4
# NorthBound     12.75

How do I crop the data using extent/crop? Once the region of interest has been cropped I would like to export it to an asc (ESRI ASCII) using the code

#writeRaster(r, filename="westAfrica.asc", format = "ascii", datatype='INT4S', overwrite=TRUE)


texts: Selecting only values that occur more than once in a field

texts:
I have a field that contains different values, some which occur more than once in the field. How do I select only the values that occur more than once?
Example

I only want to select rows that contain 'AMK' in the Subzone field as they occur more than once. How do I go about it?


texts: Unable to copy text out of MapInfo 12 info tool

texts:
As the title says, I'm currently unable to copy and paste information out of the info tool box. Is there a way that this can be done normally or do I need to write a new tool/ is there a tool out there currently that allows me do to this?

I'm using MapInfo 12 on a Windows 7 machine (if that helps).


texts: Qgis creating point within complex shapes (with holes)

texts:
I'm looking for a tool to create a point in every polygon, that will be placed in the polygon. these are complex shapes with holes. the point must not put in those holes!

the "plugin RealCentroid"  should do that, but that does not work.

I can get points, on every weight-point, but that is not always inside of a polygon!!