texts: Why does a Landsat 8 band not load into QGIS SCP plugin?

texts:
A 2015 Landsat 8 band 2 does not load into QGIS' Semi-automatic Classification Plugin window when wanting to clip multiple raster bands.

Bands are loaded from Table of Contents (TOC) as displayed above or clicking 'Refresh list'. The TOC (below) shows band 2 but, it does not appear in the SCP list. Years 2013 and 2014 do not give this problem.

2015 bands in TOC:


texts: How to write netcdf files with multiple variables?

texts:
I'm using raster stack objects and I want to use multiple stacks and write them in a single netcdf file. In case of one stack its straight forward as the below example, but how can I have both of the stacks in the same file as different variables?

library("raster")
r <- raster(system.file("external/test.grd", package="raster"))
s1 <- stack(r, r*2)
s2 <- stack(r*3, r*4)

writeRaster(s1, file="stk.nc",varname="var1",longname="var1_name", varunit="unit", xname="lon", yname="lat" ,zname="zname",format="CDF", overwrite=T,datatype='INT2S')


texts: How to limit a database layer in QGIS to the extent of another layer's feature?

texts:
I am developing a QGIS project template for a client that uses multiple PostGIS layers for viewing/editing structures that are identified from aerial images. In each new project derived from that template, users should be able to limit features from all PostGIS layers to the surrounding of  the area of the specific project, to save from loading the whole tables into QGIS, which would get slower with ever increasing table sizes. Project areas are saved in a separete PostGIS table, so it would be easy to reference from QGIS, but the layer filter does not work with project variables or similar. My first idea was to use a project variable to tell QGIS which project area we are interested in in this specific project, then query the project_area layer for its geometry, then apply a layer filter to all other postgis layers, limiting to those features that intersect with a buffer around the project area.

What would be another solution to filter database layers for intersection with another layers geometry, without manually changing layer filters or hand-writing SQL in every project derived from the template?


texts: Replacing double quotes in filename using Python?

texts:
I am working with data in an SDE, when I create a new feature class within the SDE I get a file name which contains my user name in double quotes. I do not have administrative permissions to the SDE and therefore cannot get the file name to write without having my user name in double quotes. I was wondering if there is a way to replace the double quotation marks in the file name so I can write an output to a shapefile. When I try to output to a shapefile I get an error stating 'invalid characters'. I've tried .replace('"','') but that doesn't seem to work. The only replace function I can get to work is replacing full stops in the file name by stating .replace(".","") using double quotes.

How do I replace these troublesome double quotes?


texts:
I want to use some fields in qgis form but to make it unvisible in table.
Is there any way to do it?

Now if I make the field property unvisible - I can't use it in the form.

I'll be glad to see any hints!


texts: Using FeatureCollection in Configurable Map Viewer?

texts:
I have searched everywhere online on How to use featureCollection in the operationalLayers of CMV.

Is any sample code available?


texts: dynamically aggregate data based on user input and create KML to be displayed on map

texts:
I want to dynamically aggregate data based on user input (basically firing a sql on my database ) and create a dynamic KML file out of that and display it on google map using their java script API .

How can i achieve this .

Note :All of my data is either on private cloud or in house ,initially i thought of exporting some of the data of my use (for aggregation purpose) to fusion but got lost on how to achieve my end target .


texts: Exporting NDVI time series from random points in a region

texts:
I am wondering how to export a time series of NDVI across a 20 year time series (1990 - 2010) for each point in my region. To be explicit: There are 1,000 random points in this arbitrary region. I need to export time vs NDVI for each image in this Landsat Collection for each point.
I am getting stuck with the method of extraction, because I do not want to reduce an image stack at each point. I have also seen individual point time series chart, but I would like this equivalent for each point.
This is what I have so far:

//Mask clouds in Landsat 5 imagery.
var maskClouds = function(image) {
  var scored = ee.Algorithms.Landsat.simpleCloudScore(image);
  return image.updateMask(scored.select(['cloud']).lt(20));
};
//Mask clouds and adds quality bands to Landsat 8 images.
var addQualityBand = function(image) {
  return maskClouds(image)
    //NDVI
    .addBands(image.normalizedDifference(['B4', 'B3']))
    //time in days
    .addBands(image.metadata('system:time_start'));
};
var images= ee.ImageCollection('LANDSAT/LT05/C01/T1_TOA')
  .filterDate('1990-01-01','2011-01-01');

var collection= images.map(addQualityBand);

var randomPoints = ee.FeatureCollection.randomPoints(geometry,1000);

// Add to map
Map.centerObject(randomPoints);
Map.addLayer(randomPoints, {}, 'random points',true);

Here is a link to the full script with the geometry import: https://code.earthengine.google.com/47c441094aef35980e528e3fd9bb9b0a


texts: GEE index image segmentation and classification

texts:
I have a code for performing an index calculation on chosen images and in chosen AOI. The next step I want to do is segment it by pixel value and classify them into 3 groups. I mean a range of pixel values for each group. After that from the borders/area of the 3 groups, I want to make a feature object so I can calculate indices on them.
I don't know how to do it or the options that I have for doing it.
'''// Images for the ImageCollection:

var image1 = ee.Image('COPERNICUS/S2_SR/20181003T081751_20181003T082718_T36SXA');
var image2 = ee.Image('COPERNICUS/S2_SR/20191013T081859_20191013T083011_T36SXA');
var image3 = ee.Image('COPERNICUS/S2_SR/20201007T081759_20201007T082035_T36SXA');

// Create a collection with fromImages().
var CollectionFromImages = ee.ImageCollection.fromImages(
[ee.Image(image1), ee.Image(image2), ee.Image(image3)]);
print('CollectionFromImages: ', CollectionFromImages);

// Load raster data (COPERNICUS/S2_SR).
// Load vector data (shapefile).
var fc = ee.FeatureCollection(user_account+folder+shpfile_name);
var filter = ee.Filter.inList('Name', ['Shoresh 7']);
var filteredArea = fc.filter(filter);
var S2 = ee.ImageCollection(CollectionFromImages);

print(fc)           // CRS must be WGS84 (for the shapefile)

//Image Expression for calculating BSI----------------------------------------------
var band_1 = 'B2'
var band_2 = 'B4'
var band_3 = 'B8'
var band_4 = 'B11'

var addBSI = function(image) {
  var BSI = image.expression(
  '((B11+B4)-(B8+B2))/((B11+B4)+(B8+B2))',
  {
    'B2': image.select('B2'),
    'B4': image.select('B4'),
    'B8': image.select('B8'),
    'B11': image.select('B11'),

  }).rename('BSI');
  //NDVI index
  return image.addBands(BSI.add(ee.Image(0.0)));
  //Add 0.0 value to BSI (currently not in use)
};

// Calculate BSI for images in the Image collection
var S2_withBSI = ee.ImageCollection(S2).map(addBSI);
print(S2_withBSI)


texts: How to set and configure a value relation edittype in QGIS 2.0?

texts:
I am trying to set the edit type of a vector layer to use a value relation.
The Python commands below are processed without errors or warnings, but the when I check the properties in QGIS, only the EditType is set, but the data settings are empty.
I use:

vlayer.setEditType(12, QgsVectorLayer.ValueRelation)
vlayer.ValueRelationData( 'releatedtable20131027113855772', 'key', 'value',False, True)

I got the name of the related table by inspecting:

ml=QgsMapLayerRegistry.instance().mapLayers()
dict(ml)

My guess is that I use the wrong string for referencing the layername.
I have two questions:
How to get the layername string to be used in ValueRelationData?
Should I add further code to get it all working?


texts: Why can't I paste coordinates in "Enter Map Coordinates"?

texts: Using QGIS 2.6.1-Brighton, I am attempting to use the Georeferencer. My data is from NGS datasheets. I copied and pasted the data into a text document. In the "Enter Map Coordinates" dialogue, when I click on the spaces to enter the coordinates, and click "Paste" (having copied the datum from the text document), nothing is apparent. I do not understand why the Paste function is defective in the "Enter Map Coordinates" window. Is this an error, and is there something I can do to paste my data into the spaces?


texts: Dividing the world space to geographic sub-space

texts:
I am trying to divide the world space to geographic sub-space from "general" to "specific" locations and without using latitude and longitude information. That is, I am trying to think about the world in someway like the following:

The most "general sub-space" > ... > The most "specific sub-space"

Example

Continent > State > Region > ... > City > Town > Address

What is (are) the most used "division approach(es)"?
What is (are) the most efficient "division approach(es)"?


texts: Con statement outputs using arcpy.sa.Con doesn't ignore NoData

texts:
In ArcGIS 10.2, using the Con tool, I ran a conditional statement to return RasterA where RasterB value >0. That worked fine. Transferring the code into python, the statement reads something like:

arcpy.gp.Con_sa()

BUT, I don't want to write the output to disk yet, so I'm going to use:

outraster = arcpy.sa.Con()

To see if the output is right, I save this to disk:

ouraster.save(output location)

The resulting output does not ignore NoData, it repalces it with value 0. arcpy.gp.Con_sa() successfully ignores NoData values and transfers them to the output raster as NoData. arcpy.sa.Con() fills all the NoData values with 0 in the output. This increases processing time and is unnecessary.

Any idea why this is happening?


texts:
I am from the Computer application background and want to know that which programming language is being widely used these days in the GIS domain? And what will be it's scope in the future?
How can I utilize my coding skills and contribute to the GIS community?

Also, is Node.js or Java used in this domain?


texts: qgis osx : Openstreetmap Plugin and tags

texts:
I'm using QGIS and The Openstreetmap Plugin to manage osm files.

That plugin is very useful, but limited to a few tags that will be converted into dbf file.

I'm looking for a method allowing to add some tags or rename existing in order to manage them easily in the dbf file.

I've searched for the author of this cool plugin, but no one came out.

Thank you!


texts: Convert MODIS HDF file in sinusoidal projection into GeoTIFF using Python GDAL

texts:
I'm looking for a way to convert MODIS LST data from HDF to Tiff format, and reproject from sinusoidal to WGS84 - all using Python GDAL.
An example file is here (I only care about band 0) :
http://e4ftl01.cr.usgs.gov/MOLT/MOD11A1.005/2004.03.31/MOD11A1.A2004091.h34v10.005.2007261231833.hdf
After converting all the tiles for one date, I want to use gdalmerge to merge all them into a Tiff file for the entire globe.
I've tried various combinations of gdalwarp and gdal_translate, but I think I'm not specifying the parameters properly.
Here is something I tried.
gdal_translate -of GTiff {inputfile} {outputfile}

Together with one of the following (neither worked):
gdalwarp -of GTIFF -s_srs '+proj=sinu +lon_0=0 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs' -r cubic -t_srs '+proj=longlat +datum=WGS84 +no_defs' {inputfile} {outputfile}

gdalwarp -s_srs '+proj=sinu +wktext' -t_srs EPSG:32719 -r cubic  {inputfile} {outputfile}

I'm new to GIS, though, and have been trying to figure it out from other stackexchange examples.
For anyone familiar with ArcGIS, here is a script that achieves what I'm hoping to do, but using ArcPy. It was written for MODIS ET data. Example here: ftp://ftp.ntsg.umt.edu/pub/MODIS/NTSG_Products/MOD16/MOD16A2.105_MERRAGMAO/Y2000/D193/MOD16A2.A2000193.h08v05.105.2013121202542.hdf

I want to replicate this using GDAL.

In the end I used the solution of Loïc Dutrieux for tiles away from the dateline, and for those near the dateline, I use the clipping solution of AndreJ.


texts: Does not show shortest path in networkx (in QGIS)

texts:
I used the following commands (in qgis console) to get the shortest path with astar algorithm. But the result shows 'NetworkXNoPath' every time. Why was that? Please help me to find a solution to this problem.

>>>import networkx as nx
>>>G = nx.read_shp(str(iface.mapCanvas().currentLayer().source()))
>>>len(G.edges())
18562
>>>route = nx.shortest_path(G, G.nodes()[1], G.nodes()[10])
Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "/usr/local/lib/python2.7/dist-packages/networkx/algorithms/shortest_paths/generic.py", line 136, in shortest_path
  File "/usr/local/lib/python2.7/dist-packages/networkx/algorithms/shortest_paths/unweighted.py", line 138, in bidirectional_shortest_path
  File "/usr/local/lib/python2.7/dist-packages/networkx/algorithms/shortest_paths/unweighted.py", line 203, in _bidirectional_pred_succ
NetworkXNoPath: No path between (106.8249296, -6.2047544) and (106.7920071, -6.2245956).


texts: I'm actually working on GeoNode 2.8 for a company as an internship, one of my mission is to change the basic zoom. They'd like to have a zoom on a specific region instead of having a global view of the map. For this I need to change the zoom level and the position of it but I don't know where to find this function in all these files.


texts: Cannot select / move any labels with label movement tools

texts:
I'm working within a Spatialite database.

I have within my attribute fields Created X, Y, and Rotation columns.

Originally I could select a label using the "change label" tool.

Once the data defined field types were assigned, My label movement tools appeared but I am unable to actually select any labels and move them.

Layer editing is enabled.

How do I actually move the labels, What is wrong?? From what I can tell from videos, Clicking it and "using the tool" should be enough, but as it is they are un-selectable.


texts: Spatial Interactions - ArcGIS Online - Query based on mouse-click

texts:
I have a web map, authored and deployed as a web app.
It contains a series of points and lines.
I was wondering if it was possible to create an ability to click anywhere on the map, and have the coordinate of where the click happens be the input to a spatial query of one of the existing layers?
So for example, clicking on the map produces a pop-up telling me the closest point layer, its distance and some attributes from this layer.
I was thinking of using some Arcade expressions, but it seems that you need to have the layer already created to create the popup.


texts: Converting Ellipsoidal height to Orthometric height in Python

texts:
I would like to convert Ellipsoidal height to Orthometric height in Python. I tried below:

import pyproj

EGM2008height = pyproj.Proj("EPSG:3855")
wgs84_ll = pyproj.Proj("EPSG:4326")
print(pyproj.transform(wgs84_ll, EGM2008height, 34.68016909181368, 38.31245226053967,100))

It gives this error:

> Traceback (most recent call last):   File "bos.py", line 26, in
> <module>
>     EGM2008height = pyproj.Proj("EPSG:3855")   File "C:\Users\abt\AppData\Local\Programs\Python\Python37-32\lib\site-packages\pyproj\proj.py",
> line 151, in __init__
>     super().__init__(cstrencode(projstring.strip()))   File "pyproj\_proj.pyx", line 86, in pyproj._proj._Proj.__init__
> pyproj.exceptions.ProjError: Invalid projection b'+vunits=m
> +no_defs'.: (Internal Proj Error: proj_create: unrecognized format / unknown name)

I deleted 100 (as I use ellipsoidal height) and again same error. Any help?


texts: QGIS Virtual Layer Join SQLite/ Spatialite

texts:
I am attempting to recreate a flowmap workflow from @underdark

The first step involves joining a non-spatial edges layer to a spatial nodes layer:

SELECT StartID, DestID, Weight,
       make_line(a.geometry, b.geometry)
FROM OD_all_long
JOIN nodesData a ON OD_all_long.StartID = a.ID
JOIN nodesData b ON OD_all_long.DestID = b.ID
WHERE a.ID != b.ID

The logic of this process is simple enough and I was able to achieve the same thing using R: create a line between each O-D point using geometry from node with a common ID, and then assigning the spatial data of Weight to that line from the edges dataset.

In the above code, a and b are references to the nodes table which is logical, as it is the only location with any geometry information. However, I do not understand how a and b are actually established as "references" to the node table (either through the Virtual Layers or through QGIS).

I actually get a join when I run this at the moment, but it has no geometry. As there is no further syntax explanation for make_line in the QGIS Expressions guide, I'm a little lost.

Edges data looks like:

Nodes data looks like:


texts: getting strange "grid" artifacts when creating hillshades using gdaldem

texts:
I'm trying to create hillshades from Czech elevation data (Fundamental Base of Geographic Data of the Czech Republic (ZABAGED®) - altimetry - grid 10x10 m) - a demo files is available here: http://geoportal.cuzk.cz/UKAZKOVA_DATA/GRID10x10.zip

I combine some of their txt files using some bash scripts and then create a GeoTiff using gdal_grid. The resulting GeoTiff looks like this when imported in QGIS:

As a next step I'd like to create hillshades using the Raster->Analysis->DEM and the result looks like this:

I made sure to use the bilinear option when rewarping the raster and already tried basically all available algorithms of gdal_grid.

Not sure if this is relevant, but that's how the hill shade TIFF looks like when opend in OS X Preview:

What's the source of these artifacts and how to avoid them?


texts:
I was wondering if someone could help me figure out what's going wrong with my code. It's at the very end of a code that, once a link to a geodatabase is inserted, will print the name, directory, tables, rasters, and feature classes of the geodatabase.

My problem comes with the feature class bit. For feature class, I'm supposed to also put down the field, alias and type, but it keeps popping up with the error:

  IOError: "[u'Hospitals',u'Parcels',u'Roads']" does not exist

Am I missing a line somewhere?

fcList = arcpy.ListFeatureClasses("*")

fields = arcpy.ListFields(fcList)
for field in fields:

    print("Field:       {0}".format(field.name))
    print("Alias:       {0}".format(field.aliasName))
    print("Type:        {0}".format(field.type))


texts: Can I add GeoExt or ArcGIS JS API to Openlayers + Dojo App?

texts:
I've written a basic map app with Dojo and Openlayers (using titlebar , toggle buttons , checkbox tree and some more).
I've seen ArcGIS JS Api has nice tooltip / balloon and GeoExt has nice Legend Panel which I want to add both to my app.

ArcGIS looked a little bit complicated to add (it interrupted Dojo) so I passed it.
Can I add only GeoExt just for the legend Panel or will I need to choose either GeoExt or Dojo in my app ?


texts:
What I am trying to do is average a set of points with a separate boundary.

Each of the larger circles is a fixed distance around a sample point (typically the area sampled) with each of the smaller points being a yield point. What I have been doing in the past is creating an averaged raster grid of yield and then using the point sampling tool to give me a yield per sample point (with attributes of the sample included).

While this is a fairly good process, it has hard edges and the points do not always end up in the center of the grid.

What I would like to do is use a similar process as the point sampling tool, but with a noticeable difference. I would like to return the aspects of the sampling attributes along with the average and standard deviation of the yield points within the boundary. I am only needing them in a csv/txt file to run them through some statistics to find relationships.

Is there a process out there someone has used to solve a similar problem?


texts: Return Majority from PostGIS function ST_SummaryStatsAgg()

texts:
The mentioned function

  Returns summarystats consisting of count, sum, mean, stddev, min, max for a given raster band of a raster or raster coverage.

I have an image of land cover classes (a discrete dataset), returning the mean does not make much sense since there are only classes (essentially integer numbers).

The procedure is like this: the image is clipped (using ST_clip) to the extent of a polygon (representing multiple farm fields) and then the dominant (majority) land cover class should be determined.

Question

How to get the majority from the ST_SummaryStatsAgg() function directly?

or

How to modify ST_SummaryStatsAgg() to calculate the majority?

or

Is there another function/procedure that gives me the result that I want?


texts: Zonal statistics (and other tools) will not output all records in ArcGIS for Desktop?

texts:
Using ArcGIS 10.3 for Desktop, whenever I use zonal statistics, only about 90-95% of the records actually show up in the results. The remainders are blank though, if I navigate to that polygon, there is data present.

I've tried doing this in and out of geodatabase environments, with different shapefiles and rasters as my inputs, and still the same thing. It just seems to 'max out' in some way and not return all of the records. Often, the number of records returned on the same calculation is different. The same things happens if I use the Spatial Analyst 2 toolset, or an independent toolset I downloaded meant to calculate zonal statistics for overlapping polygons.

Anyone have any ideas?


texts: Writing standalone Python scripts using PyQGIS?

texts:
I've been learning how to use the Python Console in QGIS using references from qgisworkshop.org. I'm familiar with writing standalone scripts in ArcGIS and want to learn how to do the same with QGIS.

For example, in ArcGIS 10 a simple standalone python script would be:

import arcpy

setFolder = ''

doProcess

I understand how to do this using the python console in QGIS, but I have not managed to find an example to perform the same steps using a standalone script. I suspect I've just been unlucky with my searching though. Are there any clear examples of how to do this online?


texts: Easily convert multiple attribute columns from text to integer?

texts:
I uploaded a CSV with WKT polygon coordinates. Each polygon was displayed spatially and I re-saved the file as a shapefile.

The problem is most attribute columns are being detected as text, and they should be numeric, integer values. This causes problems when sorting the data.

There are about 100 columns that are incorrectly classified as text.

Is there a way to easily convert all of them to integer?

I am using QGIS 2.18.15.


texts:
I have a layer with photos that are linked/displayed through the path attribute. Now I need to move this path and need to update it. I could of course create a whole new file but I wonder if it could be done with the field calcualtor?

Example of my path: C:/Users/Adam/Desktop/Map1/Map2/Picture.JPG

So what I would like to do is that I can simply replace section C:/Users/Adam/Desktop with for example: D:/GIS*/Map1/Map2/Picture.JPG

Is this possible in QGIS?


texts: Fixing this error "NoneType object is not subscriptable"

texts:
Since days and days, I could not figure out this issue. I got this error
'NoneType' object is not subscriptable.
lib/python3.8/site-packages/folium/features.py

The problem is at this point:
layer=folium.GeoJson(
            edges_gdf,
            tooltip=folium.GeoJsonTooltip(fields=['oneway','lanes',
            'length','speed','name'],localize=True),
            style_function=lambda x: {'color': line_color}).add_to(m)

Here is the full function:
def retrieve_data_from_postgres(network_id, Simple,
                          set_initial_crs=4326,
                          zone_radius=15,
                          intersection_radius_percentage=0.8,
                          distance_offset_percentage=0.8,
                          line_color='orange',
                          link_side='right',

                          ):
    initial_crs = "EPSG:{}".format(set_initial_crs)
    initial_crs = CRS(initial_crs)
    CRS84 = "EPSG:4326"
    default_crs = "EPSG:3857"

    edges = Edge.objects.filter(network_id=network_id)
    edges = serialize('geojson', edges,
        fields=('id', 'param1', 'param2', 'param3', 'speed', 'length',
        'lanes','geometry', 'name', 'source', 'target','network', 'road_type'))

    edges = json.loads(edges)
    edges_gdf = gpd.GeoDataFrame.from_features(edges['features'])

    nodes = Node.objects.filter(network_id=network_id).values()
    nodes_df = pd.DataFrame(nodes)
    nodes_df['location'] = nodes_df['location'].apply(geometry.Point)
    nodes_gdf = gpd.GeoDataFrame(nodes_df, geometry=nodes_df['location'])

    edges_gdf.set_crs(initial_crs, inplace=True)
    nodes_gdf.set_crs(initial_crs, inplace=True)
    nodes_gdf.to_crs(CRS84, inplace=True)
    zone_radius = zone_radius
    intersection_radius = intersection_radius_percentage * zone_radius
    distance_offset = distance_offset_percentage * zone_radius
    tiles_layer = 'OpenStreetMap'

    edges_gdf["oneway"] = edges_gdf.apply(lambda x: not edges_gdf[
        (edges_gdf["source"] == x["target"]) & (edges_gdf["target"] == x["source"])
        ].empty, axis=1)

    edges_gdf.to_crs(crs=default_crs, inplace=True)

    col_list = ['geometry', 'oneway']
    edges_gdf['_offset_geometry_'] = edges_gdf[col_list].apply(
        lambda x: x['geometry'].parallel_offset(distance_offset, link_side) if not x['oneway']
        else x['geometry'].parallel_offset(0, link_side),axis=1)

    edges_gdf.drop('geometry', axis=1, inplace=True)
    edges_gdf.set_geometry('_offset_geometry_', inplace=True)
    edges_gdf.set_crs(default_crs, inplace=True)
    edges_gdf.to_crs(crs=CRS84, inplace=True)

    latitudes = list(nodes_gdf['geometry'].y)
    longitudes = list(nodes_gdf['geometry'].x)

    # Initialize the map
    m = folium.Map(location=[latitudes[0], longitudes[0]],
                   max_zoom=18, prefer_canvas=True, tiles=tiles_layer)

    layer=folium.GeoJson(
        edges_gdf,
        tooltip=folium.GeoJsonTooltip(fields=['oneway','lanes',
        'length','speed','name'],localize=True),
        style_function=lambda x: {'color': line_color}).add_to(m)
    m.fit_bounds(layer.get_bounds())

    # Adding the nodes points
    for lat, long in list(zip(latitudes, longitudes)):
        folium.Circle(location=[lat, long], color='cyan', fill=True,
            fill_opacity=1,radius=intersection_radius).add_to(m)

    m.save('templates/visualization.html')


texts: Does QGIS use data source-specific operations?

texts:
This may be a very basic question, but I couldn't find an appropriate answer.
Does QGIS use data source specific operations where available and appropriate? For example, when using the Join attributes by location dialogue to join two layers within the same PostGIS database, will QGIS create an appropriate SQL query and run that query on the PostGIS server? Or will QGIS fetch all relevant features and do the join locally?
I know that I can write my own queries through plugins and/or DB manager, just interested in the default behaviour?


texts: How to remove features out of the extent?

texts:
I am using openlayers tile strategy and custom load function.
this.source = new Vector({
  format: new EsriJSON(),
  strategy: tile(createXYZ({tileSize: 512})),
  loader: (extent: any, resolution: any, projection: any) => {

        let source = layer.getSource();

        this.http.get(url, this.httpOptions).subscribe(
            (response: any) => {

                source.addFeatures(
                    source.getFormat().readFeatures(response)
                );
            },
            (error: any) => {
                console.error(error);
            }
        );
    }
});

But the features are icreasing in source while zooming in or out. So I want to remove features that out of the current extent from source after add. How can I do it?


texts: Point in doughnut hole incorrectly identified as within KML polygon

texts:
I have read in a KML file of CA firezones and for a certain area (around LOMPOC, CA), it seems to include anything inside a non-fire zone doughnut hole. Most of the processing is fine, just there are a few areas surrounded by fire zone where areas outside of the zones are included as being in the fire zone areas.

Is there some way to fix the issue or use another approach?

The fire threat maps are here:  https://www.cpuc.ca.gov/FireThreatMaps/

I am using the KMZ files unzipped to KML.

Code follows

import geopandas as gpd
import matplotlib.pyplot as plt
gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'

fp = "CPUC_Fire-Threat_Map.kml"
polys = gpd.read_file(fp, driver='KML')

Tier2 = polys.loc[polys['Name']=='Tier 2']
Tier2.reset_index(drop=True, inplace=True)

#test points - visually compared on map
POINT_DOUGHNUT = Point(-120.441806,34.652998)
POINT_INSIDE = Point(-120.142221,34.66407)
POINT_OUTSIDE = Point(-119.471006,36.414549)

POINT_INSIDE.within(Tier2.loc[0, 'geometry'])
> True  (correct)

POINT_OUTSIDE.within(Tier2.loc[0, 'geometry'])
> False (correct)

POINT_DOUGHNUT.within(Tier2.loc[0, 'geometry'])
> True  (wrong -- should be false)

To add the area that seems to be the problem:


texts: Inserting coordinates from Google Maps into a Database

texts:
I am using Android Google Maps to take the coordinates off a location on the map and I want to insert it into a PostgreSQL database.

I know I have to use PostGIS but how can I do that?

Who can point me in the right direction?

My insertion query in Java looks like:

String query = "INSERT INTO modul (denumire, adresa, tip_retea, geom, poza)";
query += " VALUES ('"
     + getset.getDenumire() + "', '"
     + getset.getAdresa() + "', '"
     + getset.getTipretea()+ "',
     ST_PointFromText('point(" + point.getX() + " " + point.getY()+ ")', 1), '"
     + getset.getImagine() + "');";

Which generates me the following query:

INSERT INTO modul (denumire, adresa, tip_retea, geom, poza) VALUES ('', '', '', ST_PointFromText('point(23.78977 44.320948)', 1), 'irrelevant b64 encoded image');

But how can I properly insert the coordinates into a database?

What datatype does the column 'geom' need to be?


texts: What projection is this? Is it van der Griten?

texts:
I am trying to convert some maps that came in .dra format from MapMaker to shapefiles. My problem is that the original files have no projection information, which complicates their use after exporting them using MapMaker. The map of the original data looks like this:

...which looks very similar to the van der Grinten projection, as can be seen here:

I also compared it to Robinson but that one does not look that similar. Neither does Eckert II.
One major complication is that the units of the original map do not seem to be meters, but some multiple. Specifically, if I import the world map from natural earth and compare bounds with the bounds of the exported shapefile I get different multiples, so it is not clear which unit it really is. Here's some python code that shows the issue
import geopandas as gp

# Import World shapefile from Natural Earth in WGS84
world_wgs84 = gp.read_file('https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_0_countries.zip')
world_wgs84 = world_wgs84.loc[world_wgs84.NAME.str.find('Antar')==-1].reset_index(drop=True)
world_wgs84.plot()

# Import MapMaker shapefile
mapmaker_shp = gp.read_file('https://www.dropbox.com/s/cxuoqh8msn9184o/8000.zip?dl=1')
mapmaker_shp.plot()

# Project to van der Grinten for comparison
world_vdg = world_wgs84.to_crs('ESRI:54029')
world_vdg.plot()

# compare bounds
wb = world_vdg.total_bounds # array([-19917658.38747369,  -8308978.01939191,  19974342.14947391, 14122774.74128915])
mb = mapmaker_shp.total_bounds # array([-9053.47050001, -4725.776     , 11509.66949999,  7034.509     ])
wb / mb # array([2200.00257221, 1758.2251083 , 1735.44011403, 2007.64186118])

If I define the van der Grinten projection and use different units, e.g., 1600 meters, the maps are more similar and they overlap quite closely, but not exactly as if the units in the x and y direction are different.
# Define Van der Grinten
vanDerGrinten = 'PROJCS["World_Van_der_Grinten_I", GEOGCS["GCS_WGS_1984", DATUM["WGS_1984", SPHEROID["WGS_1984",6378137,298.257223563]], PRIMEM["Greenwich",0], UNIT["Degree",0.017453292519943295]], PROJECTION["VanDerGrinten"], PARAMETER["False_Easting",0], PARAMETER["False_Northing",0], PARAMETER["Central_Meridian",0], UNIT["Meter",1600], AUTHORITY["EPSG","54029"]]'

world_vdg1600 = world_wgs84.to_crs(vanDerGrinten)
ax = world_vdg1600.plot(); mapmaker_shp.plot(ax=ax, color='r')

as shown here:
.
So the questions are

whether this is really the van der Griten projection and if so how can I transform this into WGS84 correctly
if it is not van der Griten then which projection is it?


texts:
I'm doing changes to a new field I created on the attribute table of a grid shapefile (.gpkg), the grid is circle clipped to the extent of a 4km radius buffer, it contains 499358 cells. The new field is  a short integer type and I'm assigning a number to each cell of the grid according to different types of land cover, for example, if I overlay this grid to a remote sensing image and select some cells that are 100% Forest land cover, I would put 1 in that new field, and if it is infrastructure I would put 2.

I edited the table, saved the edits and exported the grid as a .csv file but when I opened it, the new field column was completely empty and when I checked the attibute table directly on QGIS all the changes that I've donde are gone and the values are "NULL". Does anyone know what can I do about this? I'm using QGIS 3.10 long term release A Coruña


texts:
I use a custom made font for symbols on my map in ArcMap 10. Is there anyway to embed the font in the .mxd file so I don't have to distribute the font files with instructions on how to install it to my users?
On the ESRI blog I found this:

Sharing ArcGIS documents with character marker symbols
A critical thing to know is that if you want to share your ArcMap, ArcScene, or ArcGlobe document with others, they also need to have the font on their machine in order for the symbol to appear correctly.  If they don’t, the marker symbol will undergo a font substitution and will be replaced with a character from a different font. Don’t forget that you may need to check the permission on the font you want to share.

But that was from 2012, so perhaps the feature has been added since then?
Please note that I am not looking for instructions on how to embed the font into an exported .pdf - I want to keep it in the .mxd file.


texts: Using checkbox on ModelBuilder tool to control adding to data frame?

texts:
I have a geodatabase feature class and every day I want to automate the loading of my feature class into a data frame of my map.

At present in my model (below) I add into my data frame the geodatabase feature class
but I want to add a checkbox option.

i.e. choose whether to add or not the feature class using a checkbox

This is only part of the tool that I want to do.

This feature class contains reference points for an analysis of service area, which is the reason why I want to automate the data load.

Can I use a checkbox on my tool dialog to control whether, when the model runs, its output feature class gets added to the display (data frame)?


texts: Get a bearing from northing and easting data

texts:
I have a grid of points created from NOAA oceanic current data which has the currents as east and north directions that I'm needing formatted as an overall speed(got that part) as well as the direction the water is flowing in azimuth degrees.

How do I get from point A to point B?  Watched a few videos that kinda attempt to describe it but between poor English or they're just typing and music overlaid I'm not getting very far and my math really is a bit on the pathetic side.

I'm thinking it's something to do with Cosine or Tangent and I maybe need to add 90degrees based on whether the northing or easting is a negative value or not.

Currently everything's in ArcMap if there's a tool to do this...or I have the field calculator.


texts: How to delete points from CSV point layer?

texts:
I have a point layer produced from importing a Delimited Text Layer and I would like to delete some points. It seems that it is not possible to toggle editing in a point layer. How can I delete the selected points?

(QGIS version 1.8.0-Lisboa)


texts:
I have a file gdb feature class which contains cropland fields, and within the table are data fields for Landowner, Acres, and many other items.  I have a report based on the table which has no data within the detail section, but is grouped on the Landowner.  I have the summation of the Acres for each Landowner within the Landowner group footer , and a grand total of the Acres in the report footer.

What I would like is also a count of the Landowners within the report footer. (This would also be a count of the times the group footer was generated.)

Any ideas?


texts:
Is there a service that can take an address and return a polygon of the property lines (or the coordinates of the corners of the smallest possible bounding rectangle)?

The area is Oakland, CA, USA.


texts: Convert file with coordinates, keep record ID


texts: Add X Y Data not populating attribute table

texts:
I am adding X Y data to ArcMap using .csv files.

These files contain the X and Y co-ords along with a text column describing each point.

The points are added to the map correctly as an events layer, but when I view the attribute table, only the column headers are shown (see attached photo).

Can anyone suggest where I might be going wrong?


texts: Is there a reason why the quality of KMZ images, when added to a website using the Google Earth Plugin, drop in quality (it appears blurry) compared to when the same KMZ is seen through Google Earth Desktop?


texts: r - Create grid (or polygons) over shapefile that vary in area but equal in number of points encompassed

texts:
I am using R's various spatial analysis features. I have data from a number of sampling locations throughout Africa, distributed unevenly. For several reasons, it is of interest to break our Africa map down into grids and compare characteristics of grids, matched on several characteristics. Ultimately, we want to be able to randomly sample grids. However, due to the uneven distribution of sampling locations, having evenly sized grid cells is problematic, because we would end up with some grids that have only one data point and other grids that have hundreds.

Basically, I know how to generate grid cells that are equal in area, with varying numbers of sampling locations ... however, what I want is essentially the opposite. Grid cells that vary in AREA, but with equal (or, realistically, near-equal, within a certain range) numbers of sampling locations. I recognize, too, that the solution may very likely involve polygons, rather than rectangular grid cells, but I am not sure how to implement this.

Simplified, reproducible example below, which takes a shapefile of Africa, randomly generates locations across Africa (ignoring for now that many of these are randomly generated in the ocean), and layers an even grid across it:

library(maptools)
library(lattice)

# Import shapefile of Africa
data(wrld_simpl, package="maptools")
africa <- wrld_simpl[wrld_simpl$REGION==2,]

# Randomly generate 200 coordinates
locations <- cbind.data.frame( x = runif(200, min(coordinates(africa)[,1]), max(coordinates(africa)[,1])),
                               y = runif(200, min(coordinates(africa)[,2]), max(coordinates(africa)[,2])))
coordinates(locations) <- ~ x + y
proj4string(locations) <- proj4string(africa)

# Verify with plot
plot(africa)
points(locations, pch=16, bg="black")

# Create spatial grid
bb <- bbox(africa)
cs <- c(10, 10)  # cell size
cc <- bb[, 1] + (cs/2)  # cell offset
cd <- ceiling(diff(t(bb))/cs)  # number of cells per direction
grd <- GridTopology(cellcentre.offset=cc, cellsize=cs, cells.dim=cd)
sp_grd <- SpatialGridDataFrame(grd,
                               data=data.frame(id=1:prod(cd)),
                               proj4string=CRS(proj4string(africa)))

# Plot grid with randomly generated coordinates
spplot(sp_grd, "id", colorkey=FALSE,
       panel = function(..., col.regions) {
         panel.gridplot(..., border="black",col.regions="white")
         sp.polygons(africa)
         sp.points(locations, cex=1.2, pch=16, col="black")
         panel.text(...,col="red")
       })

# Can use 'over' funciton to find grid cells associated with each of our random points
# over(locations,sp_grd)

So, I know how to overlay this grid, count the number of sampling locations in each cell, extract data values for each cell, etc. etc. However, I am not sure how to proceed from here.

How would I generate grid cells (or polygons) that encompass the extent of my map area of interest, whereby the geographical area of the cells (or polygons) is variable but the number of sampling locations within each is fixed (or within a certain range)?

I am very new to the spatial analysis world.


texts: Requesting multiple features on different layers

texts:
The case: A lot of features need to be rendered on an openlayers map. Drawing all features on the map (WFS) causes performance issues (especially in IE). Therefore I have chosen to use WMS so that the client is spared the heavy rendering load.

The problem: To get information about certain features I use "Openlayers.Control.GetFeature". This works, but only when you use one layer.

I have multiple layers with different feature types. I still need to be able to use the box selection method. It may be that this exceeds the openlayers limits.

Adding multiple select controls did not seem to have any effect.


texts: dismo::randomPoints generates NA points with raster stack mask

texts:
I am using the randomPoints function in the dismo package in R to generate random points from a raster stack. The documentation for the randomPoints function says:

"The points are sampled (without replacement) from the cells that are not 'NA' in raster 'mask'."

I successfully sample random points from my raster stack rs, which consists of 9 raster layers, and put the points into a data frame for easier handling.

bg <- as.data.frame(randomPoints(mask = rs, n = 5000))
length(bg[,1])
[1] 5000

All of the points generated are valid points:

sum(is.na(bg))
[1] 0

To check whether any of these points are NA in the mask raster stack, I extract the raster stack over the points:

ext <- extract(rs, bg)
dim(ext)
[1] 5000   9

Now I check whether there are any NA's. There shouldn't be, because the randomPoints function should exclude cells of the mask raster that are NA. But still:

sum(is.na(ext))
[1] 2303

When I look through the ext data frame, I see that the NA's seem to be pretty randomly scattered through the raster stack. It's not one raster layer in particular that has all the NA's, and points that are NA in one layer aren't necessarily NA in the others. I understand that perhaps the non-NA values in my raster layers don't perfectly match up, but regardless, these should have been excluded by the randomPoints function, no?

Just to make sure, I tried this again, but this time I selected the background points using only one layer of the raster stack as the mask.

bg1 <- as.data.frame(randomPoints(mask = rs[[1]], n = 5000))

Now, when I check for NA values in the extractions, there are indeed no NA values for the raster that was used as the mask, but there are still NA values in some of the other rasters for some of the points.

It seems that passing a raster stack to the randomPoints function, instead of a single raster layer, does not result in the function excluding cells that are NA in any of the layers in the stack.

Does anyone know the best way to deal with this, if I need to have background points that are not NA in any of the raster layers in the stack? It's easily possible that there is slight misalignment in my raster layers (perhaps at land-water boundaries), but they all cover the same area and are, for the most part, well-aligned. All I can think of is to select background points, extract them over the rasters, and then exclude those that have at least one NA value, but a.) this is tedious and b.) I need to ensure that I end up with a given number of random points; since the number of NA's resulting is unpredictable, I'm not sure how to make this work.


texts: How to force ogr to treat JSON mapserver query as "Multiline" instead of "Line" geometry type?

texts:
I have a script that obtains feature data from a JSON output of an ArcGIS REST MapServer.  At present my script treats the JSON data as a "LineString" Geometry type, but, the JSON data is in fact a "MultiLineString" Geometry type.  How do I force my script to treat the data as "MultiLineString"?  My current script is:

import ogr

url = 'http://maps.six.nsw.gov.au/arcgis/rest/services/sixmaps/LPIMap/MapServer/16/query?geometry=16501408.594,-4372525.989,16523000,-4320000&geometryType=esriGeometryEnvelope&f=pjson&outFields=functionhierarchy,OBJECTID,shape,roadontype,roadnamestring,shape_Length'

ds = ogr.Open(url)
lyr = ds.GetLayer()

if (lyr.GetGeomType() == ogr.wkbLineString):
    print ('line')
elif (lyr.GetGeomType() == ogr.wkbMultiLineString):
    print ('multiline')
else:
    print ('neithor line nor multiline')

del ds, lyr

The JSON data is a MultiLineString type as it consists of seperate lists of many cordinates within a large single list.  For example, the JSON output below shows the end of one list of coordinates and the beginning of a subsequent list of coordinates and is a component of the first feature reported in the URL in the script above.

 [
   16504625.024250001,
   -4372770.7781249993
  ],
  [
   16504645.846749999,
   -4372777.386500001
  ],
  [
   16504664.026874997,
   -4372785.385625001
  ]
 ],
 [
  [
   16504684.332374997,
   -4372794.3222500011
  ],
  [
   16504706.535875,
   -4372804.0945000015
  ],
  [
   16504730.744125001,
   -4372813.5017499998
  ],

Because the script treats the data as a "LineString" different sections of the the same paths are connected by relatively long straight lines.  An example of the these erroneous relatively long straight lines are shown in the picture below which is a plot of the features in the layer as interpreted by the script.


texts: Unexpected Moran's I value in R using a chessboard example

texts:
I'm running into unexpected values for Moran's I when using the Moran.I function in the R library "ape".

My understanding is that a chessboard is an example of perfect negative spatial autocorrelation, which should yield a Moran's I value of -1.  But when I code a simple 8 x 8 grid, and give it a chessboard pattern, I get a value which, while negative, is not close to -1.

Consider this code:

library(ape)

# first, build an 8 x 8 grid
# lats and lngs will go from 0 to 7

lng <- rep(seq(0, 7, by=1), 8)

counter = 1
subCounter = 0
startNum = 0
lat = NULL

while (counter < 65) {

    if (subCounter == 8) {

        startNum = startNum + 1
        subCounter = 0

    }

    lat = c(lat, startNum)

    subCounter = subCounter + 1
    counter = counter + 1

}

# now add the black/ white chessboard pattern

chess <- rep(c(0,1,0,1,0,1,0,1,1,0,1,0,1,0,1,0), 4)

gridNum <- seq(1:64)

testGrid <- data.frame(gridNum, lat, lng, chess)

# simple Euclidean distance measure

testGrid.dists <- as.matrix(dist(cbind(testGrid$lng, testGrid$lat)))

# invert the distances, then replace the infinities with 0

testGrid.dists.inv <- 1 / testGrid.dists

diag(testGrid.dists.inv) <- 0

# now test the chessboard pattern

Moran.I(testGrid$chess, testGrid.dists.inv)

This yields the following values:

  $observed
  [1] -0.07775248

  $expected
  [1] -0.01587302

  $sd
  [1] 0.01499786

  $p.value
  [1] 3.693094e-05

The $observed value, while negative, is not the -1 I was anticipating.  Can anyone help?


texts: Get list of published items on arcgis server

texts: Is there a way to query the list of all features service layers and map layers from ArcGIS Server via ArcGIS javascript?


texts:
I am new to ArcGIS Community, and I am currently following a MOOC to learn about ArcGIS.

When I try to add Data from ArcGIS Online on the ArcMap Desktop panel, some of this data simply never show up.
It is the case for the "Terrain" map, looked for with the key word "elevation" : when I select it and click add, the map remains unchanged, so do my layers.

I tried to reboot the normal.mxt feature as suggested by someone on a similar case, but it has had little effect so far.

My version is ArcGIS Desktop 10.7.11595


texts: Handling large shapefiles and many features( e.g > 1 million points) in Mapbox and GeoServer

texts:
I have some issues. But I want to get clarity/assistance on a fundamental one.
I use GeoServer as my tile server to my client using Mapbox GL JS. Now, I have uploaded a shapefile that constitutes almost 1 million points to GeoServer. The points cover a country, like Mexico. Now when I serve this layer through mapbox, and toggle the visibility for viewing, I experience performance issues. Most times, it does not load, and causes the server to stop responding. I know there must be a way to load the tiles that would not impede performance, but I cannot figure out how.
I am concerned because in the future I might be dealing with even bigger files.
TL:DR How do I handle performance issues in Mapbox and GeoServer when I am serving shapefiles, for a region like Africa.


texts: How to solve postgreSQL error when using QGIS Cloud plugin?

texts:
I created a qgis cloud account and installed the plugin. When I attempted to Upload Data using the plugin, I got the following error:

postgreSQL exception:

ERROR:  "MD5" authentication with pgpool failed for user "database_name"

Any ideas on how to fix this error so I can upload data? Thank you!


texts: Arcgis Service Find the Closest Road Segment By giving an GPS Point

texts:
Within ArcMap, I have a road sement layer. I can add a point feature, and get the closest road segment by using near_analysis tool.

I want to do it through an ArcGIS service. Based on a given point(lan, lon), I can get the closest road segment through an Restful API.

I publish my model, my road segment and the static point feature to my local ArcGIS server. Also, I can get the closest road segment id through an Restful API, but I cannot change the point feature dynamically.


texts: Trying to use gdal2xyz in QGIS: "Error:Wrong command"

texts:
I am using the following command

gdal2xyz.dy -csv /Users/me/Downloads/esp_gpwv3_pcount_ascii_25/espp00ag.bil    /Users/documents/tttest.csv

in the commander, but I get the following error message and nothing happens "Error:Wrong command". I have installed GdalTools plugin.

Any idea why?


texts:
Simple QGIS 1.8.0 not printing to PDF. Suddenly stopped. Is this an error or something. the resultant PDFs are blank pages.

When i file print and use a pdf printer it creates blank docs. when i export to pdf it crashes. when i export to jpeg it does not print map.


texts: ArcGIS 10.0 get sde information out of .mxd

texts:
The overall goal here is to iterate through several directories of .mxd files and document which SQL Server instance of SDE they are pointing to, with the database name and feature class names.

All of our connections are direct connect, no SDE services running.

I'm trying to get the server/database info out of an mxd using ESRI's example here http://support.esri.com/cn/knowledgebase/techarticles/detail/41112, but it keeps trying to look at the .sde file that was used to create the mxd instead of getting the info straight from the .mxd.

specifically, this line

desc = arcpy.Describe(sub_str)

throws an error that it can't find the .sde file.  It's not going to since we have many users all with their own .sde connection files. The sub_str ends up being C:\Users\username\etc... depending on who created the .mxd.

Is there another way to get the datasource info out of an .mxd that doesn't rely on the .sde file? (ArcGIS 10.0)

EDIT: SOLUTION USED

On the off chance my horrible python skills can be of use to someone, this is the solution I ended up with.  I know it can be cleaned up, but I just needed quick and dirty feel free to critique:

def mxdInfo(inMxd, inCsv):

mxd = arcpy.mapping.MapDocument(inMxd)

for lyr in arcpy.mapping.ListLayers(mxd):
    if lyr.supports("SERVICEPROPERTIES"):
        servProp = lyr.serviceProperties
        lName = lyr.name
        dbName = servProp.get('Database', 'N/A')
        serviceName = servProp.get('Service', 'N/A')
        serverName = serviceName.split(':')
        userName = servProp.get('UserName', 'N/A')
    else:
        print "SERVICEPROPERTIES not supported on " + lyr.name
        continue
    if lyr.supports("DATASOURCE"):
        dSource = str(lyr.dataSource)
        fcName = dSource.split('.sde\\')
        print "Layer Name: " + lyr.name
        print "Feature Class: " + fcName[1]
    else:
        print "DATASOURCE not supported on " + lyr.name
        continue
    #write out to csv
    if not os.path.isfile(inCsv):
        csvFile = open(inCsv, 'wb')
        try:
            writer = csv.writer(csvFile)
            writer.writerow(('MXD', 'Layer Name', 'Server', 'Database', 'Feature Class', 'User Name'))
            writer.writerow((inMxd, lName, serverName[2], dbName, fcName[1], userName))
        except:
            print "error writing first row of csv"
    else:
        csvFile = open(inCsv, 'ab')
        try:
            writer = csv.writer(csvFile)
            writer.writerow((inMxd, lName, serverName[2], dbName, fcName[1], userName))
        except Exception as e:
            print "error writing to csv"
            print e
del mxd


texts: Changing SVG Search Paths without GUI at Debian | QGIS 3.10.12

texts:
I tried to load a QGIS Project which was created on windows. In this project i used a svg symbol "symbol/blue-marker.svg". In the .qgis file the path is also "symbol/blue-marker.svg".
When I start my qgisserver as request the map there are only bis ? questionsmarks because qgis can't find the blue-marker.svg.
I tried to copy the symbol folder with the symbols to /usr/lib/cgi-bin/ folder where the qgis_mapserv.fcgi file is located. Didn't work.
Then I tried to copy the symbol folder to the .qgis project file. Didn't work.
I researched a little bit and found a command to show me what the default svg paths are:
python3
from qgis.core import QgsApplication
QgsApplication.showSettings()
'Application state:\nQGIS_PREFIX_PATH env var:\t\t/usr/share/qgis/lib/qgis/server\nPrefix:\t\t\nPlugin Path:\t\t\nPackage Data Path:\t/usr/share/qgis/lib/qgis/server/share/qgis\nActive Theme Name:\t\nActive Theme Path:\t/usr/share/qgis/lib/qgis/server/share/qgis/resources/themes//icons/\nDefault Theme Path:\t:/images/themes/default/\nSVG Search Paths:\t\nUser DB Path:\t/usr/share/qgis/lib/qgis/server/share/qgis/resources/qgis.db\nAuth DB Path:\tqgis-auth.db\n'

SVG Search Paths:\t\n
From there I know that my SVG Search Paths are empty.
I can change the variable temporary with:
from qgis.core import QgsApplication
QgsApplication.setDefaultSvgPaths(['/usr/share/qgis/svg'])

but this doesn't help.
I also tried to execute this command when creating the docker image but the SVG Search Paths are also empty after that.
So in my opinion some where this variable can be set but I can't find it, also not in the QGIS Environment variables.
Here is my qgis apache2 config:
<VirtualHost *:80>
  ServerAdmin webmaster@localhost
  ServerName qgis.demo

  DocumentRoot /var/www/html

  # Apache logs (different than QGIS Server log)
  ErrorLog ${APACHE_LOG_DIR}/qgis.demo.error.log
  CustomLog ${APACHE_LOG_DIR}/qgis.demo.access.log combined

  # Longer timeout for WPS... default = 40
  FcgidIOTimeout 120

  FcgidInitialEnv LC_ALL "en_US.UTF-8"
  FcgidInitialEnv PYTHONIOENCODING UTF-8
  FcgidInitialEnv LANG "en_US.UTF-8"

  # QGIS log (different from apache logs) see https://docs.qgis.org/testing/en/docs/user_manual/working_with_ogc/ogc_server_support.html#qgis-server-logging
  FcgidInitialEnv QGIS_SERVER_LOG_FILE /var/log/qgis/qgisserver.log
  FcgidInitialEnv QGIS_SERVER_LOG_LEVEL 0
  #FcgidInitialEnv QGIS_PREFIX_PATH /usr/share/qgis/lib/qgis/server
  FcgidInitialEnv QGIS_PREFIX_PATH /usr
  FcgidInitialEnv QGIS_PLUGINPATH /opt/qgis-server/plugins
  FcgidInitialEnv QGIS_DEBUG 0

  # default QGIS project
  SetEnv QGIS_PROJECT_FILE /home/qgis/projects/world.qgs

  # QGIS_AUTH_DB_DIR_PATH must lead to a directory writeable by the Server's FCGI process user
  FcgidInitialEnv QGIS_AUTH_DB_DIR_PATH "/home/qgis/qgisserverdb/"
  FcgidInitialEnv QGIS_AUTH_PASSWORD_FILE "/home/qgis/qgisserverdb/qgis-auth.db"

  # See https://docs.qgis.org/testing/en/docs/user_manual/working_with_vector/supported_data.html#pg-service-file
  SetEnv PGSERVICEFILE /home/qgis/.pg_service.conf
  FcgidInitialEnv PGPASSFILE "/home/qgis/.pgpass"

  # Tell QGIS Server instances to use a specific display number
  # FcgidInitialEnv DISPLAY ":99"

  # if qgis-server is installed from packages in debian based distros this is usually /usr/lib/cgi-bin/
  # run "locate qgis_mapserv.fcgi" if you don't know where qgis_mapserv.fcgi is
  ScriptAlias /cgi-bin/ /usr/lib/cgi-bin/
  <Directory "/usr/lib/cgi-bin/">
    AllowOverride None
    Options +ExecCGI -MultiViews -SymLinksIfOwnerMatch
    Order allow,deny
    Allow from all
    Require all granted
  </Directory>

 <IfModule mod_fcgid.c>
 FcgidMaxRequestLen 26214400
 FcgidConnectTimeout 60
 </IfModule>

</VirtualHost>


texts: pgsql2shp ERROR: incompatible mixed geometry types in table

texts:
I have multi-linestring in .shp files and I split the lines with postgis, but that result is only geom column.

It does not works with shape file exporter

  error: incompatible mixed geometry types in table

I guess it needs at least gid, shape_leng, polyline_n ...etc columns?

How can I get data able to export .shp file with only geom column?


texts: Automatically adding vertices to network in pgRouting

texts:
I am trying to add extra vertices based on the intersections of my edges with a polygon to get access points to those polygons, from where I could run pg_drivingDistance.
Therefore, it would this way:

Create network from given roads;
Add extra vertices based on intersections of roads with polygons (they would be houses or parks, e.g.)

Is it possible to break the edge and add an extra vertex at the point of intersection without messing up the network?


texts: I think we should create a list of best desktop GIS FOSS case studies. First of all this list will be useful to anyone who considers desktop GIS FOSS for a project as an viable alternative.


texts: Finding the "bad value" in text file row

texts:
I have a text file that I am trying to bring into ArcMap to convert to a shapefile.  I am familiar with the process of bringing in the text file and then using the "Display XY Data" tool and exporting to a shapefile.

The problem I am running into is that only part of my data will display.  The first 11 records appear but the ones past that point do not.  When I try to view the attribute table, I get an error message that says "There was an error loading records.  The row contains a bad value [OIDTable11].  The row contains a bad value. [Field 5]".

I can't see any difference between the rows that appear and the row that apparently contains a bad value.  I was hoping that someone could help me find the "bad value" so that my table will import in its entirety.

Here is a screenshot showing my original text file and the rows that display in ArcMap:


texts:
I have 4 scanned maps in geotiff format and need to create a single map using OkMap with the 4 geotiff files.

How can I add adicional map files to a single map using OkMap?


texts:
I am trying to map a Shapefile (download) onto a Google Map using the Javascript API. The Shapefile describes the boundaries of NYC School Districts.

The file's projection was Lambert Conformal Conic, and so I tried to convert it to Google Mercator using QGIS. However, the coordinates still don't make sense. When I map the polygons, they cover all of Earth.

For example, one of the impossible points that I get is -8240484.27362, 4961069.9502.

How can I properly convert my Shapefile to an acceptable format?


texts:
Lets say I have a very big image With a polygon in it:

I want to polygonize this very big image with gdal polygonize. But the image is too big, so I dont have enough memory to polygonize that properly (I know there are some options but forget them for a while). So I split this image into 4 images (black lines are the borders of each image):

Now I polygonize these 4 images and I obtain 4 GeoJSON-files with 2 polygons each (one bigger white polygon, and one smaller blue polygon). Then I put all polygons into one GeoJSON so that I have a GeoJSON with 8 polygons with it. When I import it in any GIS software it looks like the first picture but when I want to choose the blue polygon, I have to click on all 4 parts. Because in fact the result looks like this (orange = border of polygon)

So I need a way to merge polygons so that neighbour-polygons of the same color get merged to one polygon.

My images are of course a little bit more complex but I had to simplify it to a circle to show the problem


texts:
I am trying to use pgrouting to find the shortest path in GeoServer. So I created new layer in sql view but its giving me error.

I created 3 table in database:

at_2po_4pgr = main table
at_2po_4pgr_noded = contains node network
at_2po_4pgr_noded_vertices_pgr= contains network topology [fields:]

and I wrote following sql statements in sql view in GeoServer

SELECT
 min(r.seq) AS seq,
 e.old_id AS id,
 e.name,
 sum(e.distance) AS distance,
ST_Collect(e.the_geom) AS geom
 FROM pgr_dijkstra('SELECT id,source,target,distance AS cost FROM at_2po_4pgr_noded',%source%,%target%,false) AS r,at_2po_4pgr_noded AS e
 WHERE r.at_2po_4pgr=e.id
GROUP BY e.old_id,e.name

but it shows this error:

Failed to create SQL view: ERROR: column r.at_2po_4pgr does not exist Position: 270


texts: Transfering a lot of features from geodatabases?

texts:
I'm having a problem where I need to export a large amount of content from 2 geodatabases of different nature (SDE -> File).

I have three options: Spatial Select, a regular SQL select and the full feature class/table.

What are the best practices to handle these large transfers?

I'm having issues with the following code:

    public void ExportAll(IGisEntity entity)
    {
        ITable table = _Context.GetEntityTable(entity.EntityName);
        IFeatureWorkspace workspace = _OutWorkspace as IFeatureWorkspace;
        ITable outTable = workspace.OpenTable(entity.EntityFeatureClass);

        if (table == null)
            throw new ArgumentException("A tabela de origem não existe.");

        if (outTable == null)
            throw new ArgumentException("A tabela de destino não existe.");

        ICursor selectCursor = table.Search(null, false);
        IRow selectRow = null;

        ICursor insertCursor = outTable.Insert(true);
        IRowBuffer buffer = outTable.CreateRowBuffer();

        int recordCount = table.RowCount(null);
        int count = 0;
        int step = recordCount;
        if (recordCount > 10000)
            step = Convert.ToInt32(Math.Round(recordCount * .1, 0));

        int remainder = 0;

        while ((selectRow = selectCursor.NextRow()) != null)
        {
            buffer = selectRow;
            insertCursor.InsertRow(buffer);
            count++;
            Math.DivRem(count, step, out remainder);
            if (remainder == 0)
                insertCursor.Flush();
        }
        insertCursor.Flush();

        CustomComReleaser.ReleaseComObject(insertCursor);
        CustomComReleaser.ReleaseComObject(selectCursor);
    }

    public void ExportSelect(IGisEntity entity, int sectorNo)
    {
        ITable table = _Context.GetEntityTable(entity.EntityName);
        IFeatureWorkspace outWork = _OutWorkspace as IFeatureWorkspace;
        ITable outTable = outWork.OpenTable(entity.EntityFeatureClass);

        if (table == null)
            throw new ArgumentException("A tabela de origem não existe.");

        if (outTable == null)
            throw new ArgumentException("A tabela de destino não existe.");

        IQueryFilter filter = new QueryFilterClass();
        filter.WhereClause = String.Format(@"{0}_GR_IDENTIFICADOR LIKE '1{1}%'", entity.EntityCode.Substring(1, 3), sectorNo.ToString());

        ICursor selectCursor = table.Search(filter, false);
        IRow selectRow = null;

        IRowBuffer buffer = outTable.CreateRowBuffer();
        ICursor insertCursor = outTable.Insert(true);

        int recordCount = table.RowCount(filter);
        int count = 0;
        int step = recordCount;
        if (recordCount > 10000)
            step = Convert.ToInt32(Math.Round(recordCount * .1,0));

        int remainder = 0;

        while ((selectRow = selectCursor.NextRow()) != null)
        {
            buffer = selectRow;
            insertCursor.InsertRow(buffer);
            count++;
            Math.DivRem(count, step, out remainder);
            if (remainder == 0)
                insertCursor.Flush();
        }
        insertCursor.Flush();

        CustomComReleaser.ReleaseComObject(insertCursor);
        CustomComReleaser.ReleaseComObject(selectCursor);
    }

When one operating on ArcSDE systems these methods tend to hang and gives me I/O errors.

I look into IGeodatabaseBridge interface and IGeoDbDataTransfer, but none of them gives me the option to export a cursor (product of a select, for example).

This is now a major issue: everytime I'm getting a RPC_ERROR_SERVERFAULT with this snippet.

        while ((row = select.NextRow()) != null)
        {
            buffer = row as IRowBuffer;
            insert.InsertRow(buffer);
            count++;
            if (count == 4999)
            {
                insert.Flush(); // this is where the error explodes at my face!
                count = 0;
            }
        }

Any choices?


texts: Changing the field colour/design on a QGIS 3 QT attribute form?

texts:
I have problems with changing the physical appearance of a custom attribute form (QGIS 3, Windows 7).

A custom attribute form can be designed, for use in QGIS, using QT / QT Designer. In QT Designer there are several methods for changing the appearance of the attribute form - changing colours of labels or fields, field background colour, fonts, and so on. In some circumstances these show up when the attribute form is used in QGIS, but in some circumstances the settings recorded in QT/QT Designer are completely ignored.

For example:
In QT designer, setting the 'stylesheet' attributes to "font-weight: bold;font-size: 18px" - specifically for a label - makes that label appear, when in QGIS, in 18px/bold.
However doing the same for a 'QLineEdit' widget only has an effect on the font weight, not the font size (it appears bold, but at a standard size).

How can the appearance of the field (i.e. the LineEdit widget) be changed? What is it which limits which elements of the style can be changed and which can't?
(Please note - I believe this is a question about QGIS, not one about QT, which is why I post it here).

EDIT -
Adding extra confusion: some settings for a QT widget design are applied to the form only when the  layer is not in edit mode. I have a QComboBox with the stylesheet setting as "font-weight:bold;font-size: 16px" which appears as large/bold only when the layer is not in edit mode.


texts: Strangely scaled and pixelated font and buttons

texts:
I'm completely new to the world of QGIS and GIS in general. I'm going through the official training manual, and I'm having fun and everything is working out all right.
For some reason, when I got back to it today after a week or so all the fonts are strangely pixelated. Also, the buttons (and the font, I'm not sure?) in the toolbar are slightly pixelated and seem scaled up a little bit. See the screenshot.
I did switch to a different monitor, but did not change settings. In the windows display settings scaling is set to 125% as it always was - I didn't change anything here either.
It also doesn't go back to normal when I use my old monitor - and QGIS is the only application I'm having this problem with (thus I'm asking here and not in the superuser forum).
Any help for that?
Edit:
I'm using Windows 10.
Also, I tried exporting a map, in the resulting pdf file the font looks good and not pixelated.


texts: Create a layer in PyQgis with the same CRS as another

texts:
I want to create a layer which use the same CRS as another layer I already have in the project. The idea is to avoid the CRS dialog.

I tried the things in this link:
How to define layer's CRS and avoid the CRS dialog in PyQGIS?

But it doesn't work because as soon as I create the layer the program launches the dialog.

qgis.utils.iface.addVectorLayer("shapefile_path", "Layer_name", "ogr")

#Use loaded layer as active layer to set the CRS
myLayer = qgis.utils.iface.activeLayer()
myLayer.setCrs(QgsCoordinateReferenceSystem(4326, QgsCoordinateReferenceSystem.EpsgCrsId))

I tried something like

crs = originalLayer().crs().authid()
layer2 = QgsVectorLayer("Point", "auxlayer", "memory")
layer2.setCrs(QgsCoordinateReferenceSystem(crs))

And it works but only after the dialog appears.


texts: CORS issue ArcGIS 10.1 IIS (Cross Origin Resource Sharing)

texts:
I have ArcGIS 10.1 installed with the web adapter for IIS. I've added a config file to the root of IIS7 to enable Cross Origin Resource Sharing (CORS) as per this page. I've pushed the cors enabled ArcGIS server into the list of default cors servers as per this page:

esri.config.defaults.io.corsEnabledServers.push("vmagstenone")).

However, when my JavaScript application which is hosted on vmagsten makes a request to the GIS server (vmagstenone) I get the error

"Origin http://vmagsten is not allowed by Access-Control-Allow-Origin.".

Edit: And the feature layer fails to load. I have also tried loading a dynamic layer and this too fails with the same error _557 (see image)

Edit: I should not be seeing this error as the above steps should mean that the server does support cross origin. This statement by esri that this error can be ignored does not apply in this case since this server should support cors. It does mean that the first error to the ESRI server can be ignored.

!

Edit: here is a sample response from a query on this layer showing that the response has Access-Control-Allow-Origin

!


texts: "Data management tools" does not exist under "Vector layer"?

texts:
I am using QGIS 2.18.12 and I need to calculate the area of overlap between two shapefiles. Help I've already found for this problem says to use "Data management tools" under "Vector layer" but the former does not exist in my version of QGIS.

Can someone explain what I need to do?


texts: OveDouble counting points that fall within 2 or more overlapping buffers in ArcGIS Desktop?

texts:
I want to analyze crime around certain types of stores in a city. I made buffers around the stores and I need to:

select and extract crimes that occurred within those buffers, and
join the two shapefiles (buffers and crimes) so that for every crime I also have information on the store it is close to (e.g. type of store, opening date).

When two buffers overlap and a crime falls within both of them, I need to double count that crime and assign it to both buffers (and both stores).

Is there a way to do this in ArcMap?


texts:
This morning, when I tried to open the Microsoft Access Database in MapInfo, MapInfo crashes.
I then tried to save a .tab file as Microfost Access Database format and the same happens, MapInfo crashes. I am currently using MapInfo version 12.5
Does anyone have any idea or experience the same issue lately?
It all worked fine up to this morning. I never had a similar issue before.


texts: Shifting cloud mask in Erdas/ArcMap/Py Script

texts:
I have the task to mask out cloud cover in approximately 20 landsat 8 scenes. Creating the cloud mask for the "white part" is easy, but I need help in addressing the clouds shadow.

General information about the shadow is that it occurs in a particular direction/distance from the cloud dependent on height/topography etc.

All I need to know is if there is a easy/automated way to shift my original cloud (white part) mask by a certain distance and direction to generally cover the shadowed area?


texts: I have shapefiles in a PostGIS database hosted in the GeoServer. I access my page with my data using localhost URL. However I need to give access to the page to other people. How can I do that?


texts: UpdateCursor based on results from SearchCursor

texts:
i'm kinda stuck with an UpdateCursor. This is what I'm trying to do: each cell of a grid stores multiple answers that were given in a survey. For each cell I must choose only ONE answer, the one that has the highest rates.

What I'm trying to do in my code is to select the rows that have the same value for the ID_GRID, compare the values in the PERCENT field, and update the FINAL_VALUE with the value from the VALUE field corresponding to the highest PERCENT value. In the example below, the answer 3310 has a rate of 76%  > FINAL_VALUE for the selected results must be 3310.

My code takes only the last value in the SearchCursor and updates the whole FINAL_VALUE field with it:

myTable = r"D:\...."
myFields = ['PERCENT', 'VALUE']

with arcpy.da.SearchCursor (myTable, "ID_GRID") as scursor:
    for srow in scursor:
        with arcpy.da.SearchCursor(myTable, myFields, "ID_GRID" + "= " + str(srow[0])) as sscursor:
        maxValue = None
        ID = None

        for ssrow in sscursor:
            if maxValue < ssrow[0]:
                maxValue = ssrow[0]
                ID = ssrow[1]

                with arcpy.da.UpdateCursor(myTable, "FINAL_VALUE") as ucursor:
                    for urow in ucursor:
                        urow[0] = ID
                        ucursor.updateRow(urow)


texts: Using QGIS comboManager to access attributes from PostGIS

texts:
I have been tasked with creating a plugin in QGIS that would connect to a PostGiS database, then let the user (based on its user ID) read the attributes of a table for a particular field, then choose based on this selected attribute, another attribute and then add the view to QGIS.

I created my plugin using Plugin Builder.

Then I used QgsCredentials and QgsDataSourceURI to connect to my database. Works fine. I get the users to enter their username and password to connect. It looks like this:

Once this step is passed, the user should be able, via combobox, to select the name of their harvest (I am working with gardeners) then the id of the harvest they want to update. The structure I created looks like this:

Now, I have been trying to use the QGIS ComboManager and using this post and this blog post (in spanish) as references. This option works perfectly fine when my layers are in the TOC.
However, using ComboManager, I feel like it would work only if I have a vector layer in my TOC. But the users won't have access to the full table, only to a spatial view. Maybe QGIS ComboManager is not the right approach.

I get this error:

def_init_(self, widget, comboBox, comboBox_2, initField="NOM_RECOLTE_", options={}):
SyntaxError: invalid syntax

It doesn't seem to like my initField. I checked on this link but I cannot identify my mistake. Does someone can put me on the right track?


texts: Displaying weighted counts through symbology

texts:
I have an Excel table of observational data.

Loc         Latitude    Longitude       Date       Survey?    Observed?
YDL         42.20833333 -83.62344444    8/15/16    NO         YES
YDL         42.20833333 -83.62344444    8/16/16    NO         YES
YDL         42.20833333 -83.62344444    8/17/16    YES        YES
YDL         42.20833333 -83.62344444    8/18/16    YES        NO
Courthouse  42.24027778 -83.61861111    8/15/16    YES        NO
Courthouse  42.24027778 -83.61861111    8/16/16    YES        YES
Courthouse  42.24027778 -83.61861111    8/17/16    YES        NO
Courthouse  42.24027778 -83.61861111    8/18/16    YES        YES
EMU         42.24027778 -83.61861111    8/15/16    NO         NO
EMU         42.24027778 -83.61861111    8/16/16    NO         NO
EMU         42.24027778 -83.61861111    8/17/16    NO         NO
EMU         42.24027778 -83.61861111    8/18/16    YES        YES

I have converted this to a feature class and added it to a map:

I want these points to be proportional symbols, based on the count of 'Yes' values in the 'Observed' column. For example, YDL has 3 YES values, so it should be larger than EMU (2 YES) and EMU (1 YES). I managed to do this by basing the symbology off the Observed column, and changing the symbols based on the 'Count' column in the symbology tab of each category:

Here's the catch. I want to somehow weight those counts against the number of values in the SURVEY column, because we surveyed some locations more frequently than others. For instance, we only surveyed once at EMU, but we surveyed 4 times at Courthouse, which gives a higher probability of observing something at the Courthouse.

Can I do this in ArcMap? It would be nice to not have to restructure my data,  but maybe it would be easier if I just had one row for each point, and a separate column to create a weighted count as described above?


texts: Converting points from MySql database to WKT lines using PHP?

texts:
In MySql database I have points (pairs of latitudes and longitudes represented as doubles) like this:

LATITUDE | LONGITUDE
---------|----------
45.045664| 46.084528
26.938933| 60.903848
       .....

Now, what I need is to get these points from my database (which is not a problem at all) and echo them in a WKT format or WKT line (this part I don't know). This PHP script is going to be a part of an API used by Windows App.

How can I convert points to WKT format?

Perhaps there is a plugin I can use.


texts: Creating custom proj4 string from components

texts:
I'm trying to work with NOAA's NDFD weather data, which are released as rasters in GRIB format. Specifically, I'm currently trying to get the pixels that intersect an arbitrary geometry.

Information about the grid's projection is here, and an example GRIB file is here (42MB).

When I load the file with rasterio, I see that the bounds are essentially 0-1377 in the y direction and 0-2145 in the x direction:

import rasterio
r = rasterio.open('YEUZ98_KWBN_201701010519')
r.bounds
>>> BoundingBox(left=-0.5, bottom=1376.5, right=2144.5, top=-0.5)

When I take that same file and run gdalsrsinfo on it, I get a proj4 string:

> gdalsrsinfo YEUZ98_KWBN_201701010519 -o proj4

+proj=lcc +lat_0=25 +lon_0=265 +lat_1=25 +lat_2=25 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs

However if I take an arbitrary point in the US, and try to project it to that projection, I get x and y values way outside the 0-1377, 0-2145 range I expect:

import geopandas as gpd
import pyproj
from shapely.geometry import Point

proj_str = '+proj=lcc +lat_0=25 +lon_0=265 +lat_1=25 +lat_2=25 +x_0=0 +y_0=0 +R=6371200 +units=m +no_defs'
crs = pyproj.CRS.from_proj4(proj_str)

point = Point(-122.4194, 37.7749)
gdf = gpd.GeoDataFrame(geometry=[point], crs={'init': 'epsg:4326'})
gdf = gdf.to_crs(crs)
list(gdf.iloc[0].geometry.coords)

# [(-2456716.706510244, 1682088.477346224)]

I expected this code to reproject the WGS84 point to a point on the NDFD grid, but it doesn't. Why not?

Edit: output of gdalinfo:

> gdalinfo YEUZ98_KWBN_201701010519
Driver: GRIB/GRIdded Binary (.grb, .grb2)
Files: YEUZ98_KWBN_201701010519
Size is 2145, 1377
Coordinate System is:
PROJCRS["unnamed",
    BASEGEOGCRS["Coordinate System imported from GRIB file",
        DATUM["unnamed",
            ELLIPSOID["Sphere",6371200,0,
                LENGTHUNIT["metre",1,
                    ID["EPSG",9001]]]],
        PRIMEM["Greenwich",0,
            ANGLEUNIT["degree",0.0174532925199433,
                ID["EPSG",9122]]]],
    CONVERSION["unnamed",
        METHOD["Lambert Conic Conformal (2SP)",
            ID["EPSG",9802]],
        PARAMETER["Latitude of false origin",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8821]],
        PARAMETER["Longitude of false origin",265,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8822]],
        PARAMETER["Latitude of 1st standard parallel",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8823]],
        PARAMETER["Latitude of 2nd standard parallel",25,
            ANGLEUNIT["degree",0.0174532925199433],
            ID["EPSG",8824]],
        PARAMETER["Easting at false origin",0,
            LENGTHUNIT["Metre",1],
            ID["EPSG",8826]],
        PARAMETER["Northing at false origin",0,
            LENGTHUNIT["Metre",1],
            ID["EPSG",8827]]],
    CS[Cartesian,2],
        AXIS["easting",east,
            ORDER[1],
            LENGTHUNIT["Metre",1]],
        AXIS["northing",north,
            ORDER[2],
            LENGTHUNIT["Metre",1]]]
Data axis to CRS axis mapping: 1,2
Origin = (-2764474.350731994025409,3232111.710792394354939)
Pixel Size = (2539.702999999999975,-2539.702999999999975)
Corner Coordinates:
Upper Left  (-2764474.351, 3232111.711) (130d 7'22.60"W, 49d56'49.26"N)
Lower Left  (-2764474.351, -265059.320) (121d33'48.72"W, 20d10'43.04"N)
Upper Right ( 2683188.584, 3232111.711) ( 60d51'57.95"W, 50d 6'46.52"N)
Lower Right ( 2683188.584, -265059.320) ( 69d11'54.73"W, 20d19' 6.40"N)
Center      (  -40642.883, 1483526.195) ( 95d27' 8.65"W, 38d13' 5.87"N)


texts:
I want to make a hovmoller diagram; I want to change the order of Y axis (to start with 2010 - on bottom and finish 1970 - to be on top:

library(zoo)
library(rasterVis)
library(raster)

download.file('http://iridl.ldeo.columbia.edu/SOURCES/.CAC/.sst/data.nc',
              destfile = 'SST.nc') # here you can find the data (about 3,9 MO size)

SST <- brick('data.nc')
idx <- seq(as.Date('1970-01-01'), as.Date('2003-03-01'), by='month')
tt <- as.yearmon(idx)
SST <- setZ(SST, tt)
names(SST) <- as.character(tt)

hovmoller(SST,
          panel = panel.levelplot.raster,
          interpolate = TRUE,
          yscale.components = yscale.raster.subticks,
          par.settings = BuRdTheme)

So I need to have 1970 on top of the graph instead of 2010. Anyone can offer some help please?


texts: ArcPy Data Driven Pages AttributeError: PageLayoutObject: Error in exporting pages

texts:
I am trying to create an atlas with embedded title page and extra reports pages. Seems that I need to use Python scripts for this. I am not programmer and never worked with python. I found necessary script but then faced with a problem to make an export of pages from mxd

Please see my script:

        import arcpy, os

# Create an output location variable
outDir = "K:/IM/02_GIS/01_Projects/Final"

# Create a new, empty pdf document in the specified output location folder
finalpdf_filename = outDir + "/Test_Atlas_Final.pdf"
if os.path.exists(finalpdf_filename):
  os.remove(finalpdf_filename)
finalPdf = arcpy.mapping.PDFDocumentCreate(finalpdf_filename)

# Add the title page to the pdf
finalPdf.appendPages("K:/IM/02_GIS/01_Projects/TitlePage.pdf")

# Add the overview map to the pdf
finalPdf.appendPages("K:/IM/02_GIS/01_Projects/IndexMap.pdf")

# Export the Data Driven Pages to a temporary pdf and then add it to the
# final pdf. Alternately, if your Data Driven Pages have already been
# exported, simply append that document to the final pdf.
#
mxdPath = "K:/IM/02_GIS/01_Projects/AOR_General_v2_Atlas_Test1_1.mxd"
tempMap = arcpy.mapping.MapDocument(mxdPath)
tempDDP = tempMap.dataDrivenPages
temp_filename = "K:/IM/02_GIS/01_Projects/tempDDP.pdf"
if os.path.exists(temp_filename):
  os.remove(temp_filename)
tempDDP.exportToPDF(temp_filename, "ALL")
finalPdf.appendPages(temp_filename)

# Insert the pdf pages containing the reports and graphs into the final pdf
#
finalPdf.insertPages("K:/IM/02_GIS/01_Projects/Report_pg4.pdf", 4)
# Update the properties of the final pdf
finalPdf.updateDocProperties(pdf_open_view="USE_THUMBS",
pdf_layout="SINGLE_PAGE")

# Save your result
finalPdf.saveAndClose()

# Delete variables
del finalPdf

Error I keep getting

Runtime error
Traceback (most recent call last):
  File "<string>", line 28, in <module>
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\utils.py", line 182, in fn_
    return fn(*args, **kw)
  File "c:\program files (x86)\arcgis\desktop10.4\arcpy\arcpy\_mapping.py", line 447, in exportToPDF
    return convertArcObjectToPythonObject(self._arc_object.exportToPDF(*gp_fixargs((out_pdf, page_range_type, page_range_string, multiple_files, resolution, image_quality, colorspace, compress_vectors, image_compression, picture_symbol, convert_markers, embed_fonts, layers_attributes, georef_info, jpeg_compression_quality, show_selection_symbology), True)))
AttributeError: PageLayoutObject: Error in exporting pages

It works well till this step # Add the overview map to the pdf
It adding overview map perfectly and error occurs during proceeding of the next step.


texts: Converting shapefile to KML using JavaScript or PHP?

texts:
My client needs to generate maps from shapefiles. Through my initial search, I came to know that we need to convert shapefiles to .kml then we can generate maps from google maps. So I searched for a script but could not find anywhere regarding that. So please help me with this problem. I am new to this, tell me am I going in right way?

One more thing, can we generate shapefile lines on a jpg image on the right place if I take an image of a state?


texts: How to achieve equivalent to Dissolve from ArcGIS for Desktop using ST_Union?

texts:
I am trying to achieve an equivalent to the Dissolve tool from ArcGIS for Desktop using ST_Union from PostGIS but it seems I am not getting expected result.

I have one table which has certain attributes with the Polygon Geometry. (like FID, Locstat, Loccnt, Shape)

Here is my query :

SELECT c.fid, ST_Union(c.boundaryshape) FROM c Group by c.fid,c.boundaryshape;


texts:
I used the function landgenreport from the package popgenreport in R to create an object result.

This object is a list with various objects including SpatialLines class objects. I am trying to write these lines to a file readable in GIS software such as ArcGIS or QGIS. The problem is that I can not find a way to do that. I've tried several functions from several packages and I always get an error.

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

results$leastcost$paths[[1]][[78]] #this is the specific SpatialLine
linha78<-results$leastcost$paths[[1]][[78]]
writeSpatialShape(linha78,"linha78.shp")
Error in writeSpatialShape(results$leastcost$paths[[1]][[78]],"linha78.shp"):x is aSpatialLinesobject, not a compatible Spatial*DataFrame

So I tried to transform the SpatialLines object to a SpatialDataFrame and got another error message.

x = data.frame(linha78)
Error in as.data.frame.default(x[[i]], optional = TRUE) :
cannot coerce class "structure("SpatialLines", package = "sp")" to a data.frame

I've tried a lot a different possible solutions I read here and there, but it simply does not work. As such:

writeOGR with a spatialpolygon simplified by gSimplify

How to export a dataset with "SpatialPolygonsDataFrame" as a shapefile

http://r-sig-geo.2731867.n2.nabble.com/Saving-a-spatial-point-object-into-a-shapefile-td5704759.html


texts: How to make a WFS request to GeoServer with OpenLayers 3

texts:
This is OpenLayers 2 code I'm upgrading to OpenLayers 3 (OL3).
A set of filters (for example fromDate and toDate to for query) are created and inserted. How can I do this in OpenLayers 3?

I have a set of parameters which the user changes before the query. Filter is not a part of OL3, so how can I create such request to GeoServer with my parameters? Can anyone provide an example of OL3 syntax?

       // create filter
        var dateFilter = new OL.Filter.Comparison({
        type : OL.Filter.Comparison.BETWEEN,
        property : "stdate",
        lowerBoundary : date.fromDate,
        upperBoundary : date.toDate
    });

        var filter = ns.query.constructFilterString(filterBbox, date, attr);

This is where the request is created.

    // GetFeature request with filter, callback handles result from server
       ns.WebFeatureService.getFeature({
        TYPENAME : "floats",
        FILTER : filter
    }, displayFeatures);

    // fires a GetFeature WFS request to server
       function getFeature(callback) {
    // some default parameters that will be set automatically if not
       var parameters = {
        REQUEST : "GetFeature",
        SERVICE : "WFS",
        VERSION : "1.0",
        OUTPUTFORMAT : "json",
        TYPENAME : "floats",
        FILTER : filter
    };

      asyncGetRequest(parameters, callback);
}


texts:
(I gave up. I already spend 10 hours reading internet and I am not any closer to get the job done.)
If I start at the beginning:

I have country map divided into 12 regions.
I have country map divided into 56 areas.
I have country map divided into 211 townships.
I have country map divided into 2700 local areas.
I have buildings map approx. 500 thousands records.

I want to build interactive map, but on old computers OpenLayers(geojson) with 211 records is kind of slow.
My idea was to show user only regions, then a user select a region and show areas for selected region, user select an area, and then I show user only townships for selected area and so on.
I have no problem manually connect townships to areas and areas to regions, just need to split them.

My first idea was to split polygons in layer (Areas) into separate files with QGIS and load correct file in OpenLayers.

But then when searching the web how to auto split polygons in layer and save them in separate files I asked myself should I use Geoserver or something similar.

I hope I described my problem good enough to get some advice, should I use GeoServer or use separate files and how can I save polygons in layer to separate files?
How should I approach the problem?

Thanks.


texts: Programmatically naming outputs from the QGIS Modeler?

texts:
I am attempting my first use of the modeler from the Graphical Modeler from the Processing menu in QGIS (2.18.3). I made a simple model for using meanshift segmentation. It just has a raster input and the algorithm for the segmentation.

I have several small rasters that I want to segment and so I want to run this iteratively through the rasters. I am not sure how that is done, so I am trying to run the model as a batch process. The batch parameters require me to name an output for each line. It seems tedious to have to write out a path and name for each output. I don't know how to simplify this. In ArcGIS' model builder, I add a generic name with an iterator to just number each new file. The iterator in that modeler is %i%. How should I use the process model to iterate through the list? Is batch processing the right way? How do I designate the names for the output without typing each out individually? I'd be fine if the process just referenced back to the raster names.

EDIT after Oto's answer. I don't get that menu from .... I only get a Save window as shown. You can see that I clicked in the .... I right-clicked and left-clicked and didn't get the save options that you show in your answer.


texts: I can't find the right python.exe in QGIS LTR 2.18 - OSGEO464W

texts:
I have recently installed QGIS 2.18.19 LTR with the OSGEO464W installer. My system is Windows 7 64bits. I installed Eclipse and PyDev successfully. However when configuring PyDev I have to indicate the 'Interpreter Executable' which in the tutorial I am following would be C:\ProgramFiles\QGIS2.18\bin\python.exe.

If I search the file, I find it in C:\OSGEO4W64\apps\Python36\python.exe, however I am certain that I am running Python 2.7.5  when I check in the python console from inside QGIS and not Python 3.
My question is: where can I find the executable or the file that opens the right python 2.7.5, what is the name that I have to provide for the PyDev Interpreter.

Below an screen capture


texts: ArcPy UpdateCursor auto pausing after every 1000 records

texts:
I have an update cursor running that updates a field based on a lookup table (Python dictionary).  The cursor completes fine, but seems to pause and hold after each 1000 records - I have a test in it to write out the time each 100 records, see code.  If you look at the results - you can see the first 1000 records don't seem to take any time, but then the 1100th result is delayed.  I am running this on an Oracle DB.

x = 0
cur = arcpy.da.UpdateCursor(property_fc,["FACILITYID","GLOBALID"])
for row in cur:
    x = x + 1
    if (x % 100) == 0:
        print("At {}, {}".format(x,datetime.now() - vST))
    if row[1] in global_id_lut:
        row[0] = global_id_lut[row[1]]
        cur.updateRow(row)


texts: What is the right projection for area calculation in Kenya?

texts: I want to calculate the area of counties (1st admin level) in Kenya in ArcGIS. What is the right projection for that? I know I should use an equal area projection, but when I use them (I tried six different) it varies heavily to the county areas given in wikipedia. I could also use UTM37, but by definition it's not an equal area projection.


texts: In Leaflet style function, option to prevent display completely?

texts:
In a Leaflet style function, is there an option which completely prevents the display of a feature? In other words (pseudocode):

function styleFunc(feature) {
    if (feature.properties['blah'] == 1) {
        return {
            color: 'red';
        }
    } else {
        return {
            DISPLAY: NONE;
        }
    }
}

What can I use instead of DISPLAY: NONE to prevent the given features from rendering? I don't just want to make them transparent, since I guess then popups etc will still operate.


texts: Convert mapinfo styles toward qgis styles

texts:
Up to now, I has worked with mapinfo. I have a big vector layer (points) with a complex symbology.

For several days, I gave up mainfo for qgis. Is there a way to converting my mapinfo files (.tab) toward a file format supported by qgis. of course, i would like saving the symbology (styling) during the file conversion.

I would like convert a big point layer from .Tab toward .SHP and saving the style into Qgis. The symbology of my layer is quite complexe (variables of coulours, shape, symbols). I'm working with a windows 7 32 bits. Mapinfo 12.0 and qgis 2.6 are installed.

I downloaded the plugin from this following repository :

https://github.com/NathanW2/MapInfo-to-QGIS-style-generator

I opened the prompt command OsGeo4W shell available with my qgis 2.6.

I tried to execute the plugin. I received an error message. This message seems to be dealing with a fied which must be supplied (with -c argument) but i don't understand.

Which column must be supplied?

The current structure of my command is :
python mapinfoToQgis.py pathINPUT.Tab pathOUTPUT.qml --UseMapInfo

I paste two screen shots :

A. focus on my .Tab mapinfo POINT layer with complex symbology

B. error mesage via prompt command

Could you throw light for me?


texts:
I try to deploy GeoServer cluster with GeoServer version 2.16.x.

I try to deploy my cluster in Ubuntu 16.04.

I have 2 VM machine with IP

IP1: 203.162.141.7
IP2: 203.162.141.13

Here is my install step (I perform the following steps on both machines)

Install default-jre and default-jdk
I dowload Geoserver package "Platform Independent Binary" by

wget https://nchc.dl.sourceforge.net/project/geoserver/GeoServer/2.16.0/geoserver-2.16.0-bin.zip

Then I unzip this file geoserver-2.16.0-bin.zip, put it in this path

/usr/share/geoserver/geoserver-2.16.0/

Then I download 2 file:

geoserver-2.16-SNAPSHOT-activeMQ-broker-plugin.zip
geoserver-2.16-SNAPSHOT-jms-cluster-plugin.zip

Unzip anh put it in this path

    /usr/share/geoserver/geoserver-2.16.0/webapps/geoserver/WEB-INF/lib

4. Start geoserver

root@ubuntu16:~# echo "export GEOSERVER_HOME=/usr/share/geoserver/geoserver-2.16.0" >> ~/.profile
root@ubuntu16:~# . ~/.profile
root@ubuntu16:/usr/share/geoserver/geoserver-2.16.0# cd bin/
root@ubuntu16:/usr/share/geoserver/geoserver-2.16.0/bin# sh startup.sh

Result:

Geoserver start like this:

when I started two machines simultaneously, it show that error:

18 Nov 23:37:18 WARN [network.DiscoveryNetworkConnector] - Could not start network bridge between: vm://4e2acd70-0c58-4ca5-a91e-02ae8a6730ad?async=false&network=true and: nio://localhost:36433 due to: java.net.ConnectException: Connection refused

I don't know how to fix this and my install step to deploy cluster is correctly or not?


texts: QGIS crashed when using labelling function

texts:
I am using QGIS 2.18.2 on a MAC OS Sierra.

Every time I try to "show labels for this layer" in the layer styling window, my program crashes as soon as I define the layer (by selecting "label with"). I see from a previous post that someone else had this same problem with an earlier version of QGIS, but that problem went unresolved in the comment string. Wondering if a solution has been found yet?


texts: Deleting duplicate points from two different layers

texts:
I have two layers with the same feature ("Object ID") and same location

All lightning strikes
Cloud strikes

I want to delete the 'Cloud strikes' that are still to be found in the 'All lightning strikes' based on the location or the same "Object ID". How can I do that?

Notice, these duplicates are NOT in the same vector layers, so I cannot use v.clean or Remove Duplicates as they only take one layer into account.

And, I do not want to merge the files and then delete as my QGIS messes up my layer projection.

Edit: The layers are the same projection, and the "Object ID"s are features in both attribute tables.