texts: How to handle AutoCAD map to QGIS workflow?

texts:
I'm looking for some assistance and like many here, I am pretty much a GIS novice.

I'm working on a project to map our local cemetery (about 40 acres), locating sections, lots, and graves, etc.

I was able to obtain an aerial photo of the cemetery and, using a series of about 6 reasonably accurate GPS waypoints, create a georeferenced TIFF using QGIS and GDAL.  This seems to have worked fairly well. Imported as a layer into QGIS, it overlays the Google imagery from the Open Layers plugin well enough for my current use.

Since I have a little more drawing control, I've been trying to use AutoCAD map for creating the features (call this the “master.dwg”).   I imported the TIFF file and drew a few features based on the TIFF image.  Exported them to shapefiles, opened in QGIS and this looked pretty good as well (with the previous TIFF and OpenLayers shown).

My next step was to create a layout of one of the sections in AutoCAD map.  I created a drawing using architectural units so I could enter values as 25’ 6”, etc.  Once I completed the layout, I tried inserting the section into the master.dwg.  For some reason, the scales seem to be massively out of whack.  Basically, the inserted section is probably hundreds of times larger than the georeferenced TIFF.  The large dimension of the section is only about 200 feet, yet when zoomed to extents, the TIFF file doesn't even appear.  It is necessary to zoom in quite a bit before it even appears as point in the drawing.

Using the measure tool in AutoCAD, measuring the distance between two known points on the TIFF literally registers a Zero (I would have expected it to be about 200 ft).  In QGIS, it measures about as expected.

Any thoughts as to what is happening here? I'm at a loss right now.

===================================================

I've made a little progress.  Today, I was able to get some shapefiles from our local GIS dept. of the area around the cemetery.  I've tried MAPCONNECT and MAPIMPORT with the same results...   I create a new drawing, assign the projection (MI83-SIF) and then attach / import the shapefiles. Units for the drawing are standard decimal.  Using the MEASUREGEOM, I measure between two points and get 204 (which is what I would expect in feet).  When I change the units to Architectural, then MEASUREGEOM returns just over 17'.  So, it appears that, while "204" in the shapefile should be feet, AutoCAD Map is interpreting it as inches.


texts: Does anyone know if a mobile version of QGIS is planned for the Blackberry Playbook. I know there is an android version, I tried converting it to .bar and gave up when it wouldn't load because it needed Minstro to configure. Any help here would be appreciated.


texts:
I want to be able to loop through the layers in my map, and return the names of layers that  are time-enabled.
I understand that I will need to check TimeInfo and this could be enabled on a dynamic, feature or Image service within the map.
I am looking for the tidiest way of doing this, in its own function that can be called and a name of a layer returned.  I have been messing about with some If statements within a 'for each' loop, but I keep struggling and making it more complicated than it needs to be.
Any advice/pointers/snippets on how to code this up in Flex would be appreciated.


texts:
I am using SplitLayerByAttributes toolbox. It’s a very useful tool for me to process certain activities. As you know, I can split the feature class as per 1 field. What if I want to split the features by 2 fields.
In the below mentioned example:
Here I would like to split the feature based on 2 fields. i.e “parcel_type” and “priority_takes”
Main field is “parcel_type”
Secondary field is “priority_takes” = Corridor

Right now I have run the toolbox twice with 2 different field and you can see the results over here. If at all possible, Splitting a feature with 2 fields then the corridors will be the subdivision in the “Parcel_types” features. Please let me know if you have any solution or any other toolbox available in arcscripts website.


texts: How to make a Python Multiprocessing Script work on Windows OS

texts:
I use the Python multiprocessing module to speed up a function (get_craters_for_buffering_BCC) in my script. In this function, I want to identify points (passed as a list: crater_features_list) within a certain geodesic distance from a polygon (passed as a polygon: geom) using OGR and Shapely.

The idea is to split the original list (crater_features_list) into n parts (crater_features_list_splitted[n]), to pass the respective parts to individual processes (processes[n]) and to pass the results from the function using multiprocessing Queues (out_q, out_q_2, out_q_3).

The procedure works fine when I execute the script in Linux OS. However, when I execute the script in Windows OS, I get the following error message:

Traceback (most recent call last):
  File "C:\Christian\Projects\2016-07 Geodesic Polygon Buffer\jo\Geodesic_Buffering_Areas18_2-Multicore.py", line 879, in <module>
    processes[process_count2].start()
  File "C:\Python27\ArcGIS10.4\lib\multiprocessing\process.py", line 130, in start
    self._popen = Popen(self)
  File "C:\Python27\ArcGIS10.4\lib\multiprocessing\forking.py", line 280, in __init__
    to_child.close()
IOError: [Errno 22] Invalid argument

This is my script:

def get_craters_for_buffering_BCC(out_q, out_q_2, out_q_3, crater_features_list_part, geom):

    """ ... determination of points within a certain distance from geom Polygon: The results
    are saved in craters_inside_area and craters_within_range lists ... """

    craters_for_counting_list = craters_inside_area + craters_within_range

    # passing information to main script using Queues

    out_q.put(craters_for_counting_list)
    out_q_2.put(len(craters_within_range))
    out_q_3.put(len(craters_inside_area))

if __name__ == '__main__':

    multiprocessing.freeze_support()

    # get geometry from file

    path_to_area_file = "path/to/shapefile.shp"
    driver = ogr.GetDriverByName('ESRI Shapefile')
    area_file = driver.Open(path_to_area_file)
    layer = area_file.GetLayer()
    area_feature = layer.GetFeature(1)
    geom = area_feature.GetGeometryRef()

    # crater_features_list was created at a different point of the script.
    # Here, the list is splitted in 20 parts for multiprocessing
    # (or 21 parts - depending on the number of features)

    crater_features_list_splitted = [crater_features_list[v:v + (len(crater_features_list)/20)] for v in xrange(0, len(crater_features_list), (len(crater_features_list)/20))]

    out_q = multiprocessing.Queue()
    out_q_2 = multiprocessing.Queue()
    out_q_3 = multiprocessing.Queue()
    processes = dict()
    process_count = 0

    # definition of process tasks (processes[n]) for the individual lists

    for crater_features_list_part in crater_features_list_splitted:
        processes[process_count] = multiprocessing.Process(target = get_craters_for_buffering_BCC, args=(out_q, out_q_2, out_q_3, crater_features_list_splitted[process_count], geom))
        process_count += 1

    # execution of multiprocessing processes

    process_count2 = 0
    for crater_features_list_part in crater_features_list_splitted:
        processes[process_count2].start() # This is where the error occures
        process_count2 += 1

The error occures only on Windows OS when I try to start the multiprocessing process:

processes[process_count2].start()

I have only little experience using multiprocessing. I tried different workarounds such as using multiprocessing.Pool, passing only string arguments to the function and calling the script from the Windows command line. They all resulted in the same error message.
How can I use the multiprocessing script on Windows OS? I use Python 2.7 (32 bit on Windows, 64 bit on Linux).


texts: Are there open standards for storing multi-modal routing graphs?

texts:
I'm looking for any open standards (or maybe best practices) for storing multi-modal routing graphs, eg. info for pedestrian, bike, public transport and car routing, etc.

Every data provider seems to have their own format. My goal is to convert data from different providers into a common standardized format and then build applications from there.


texts:
I am trying to perform PCA on MODIS imagery. I got my code to run using MODIS (it also worked for Landsat), but when I mask my MODIS collection the PCA fails. This had to do with missing values, as GEE doesn't run on spares matrices.
Questoin: How do I mask bad pixels during a PCA? One can not mask before the PCA, so how can I mask pixels after PCA?
Here is my UPDATED code. Extent is my ROI.
var dataset2019 = ee.ImageCollection('MODIS/006/MCD43A4').filter(ee.Filter.date('2018-12-01', '2019-04-20'));
var dataset2020 = ee.ImageCollection('MODIS/006/MCD43A4').filter(ee.Filter.date('2019-12-01', '2020-04-20'));
var modisCol = dataset2019.merge(dataset2020)

// Create a QA mask function
var filterqa = function(image){
  var mask1 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band1").eq(0);
  var mask2 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band2").eq(0);
  var mask3 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band3").eq(0);
  var mask4 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band4").eq(0);
  var mask5 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band5").eq(0);
  var mask6 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band6").eq(0);
  var mask7 = image.select("BRDF_Albedo_Band_Mandatory_Quality_Band7").eq(0);
  return image.updateMask(mask1).updateMask(mask2).updateMask(mask3).updateMask(mask4)
  .updateMask(mask5).updateMask(mask6).updateMask(mask7);
};

// Prepairing imagery for PCA
var PCA = function(image){
  var image = image.select(['Nadir_Reflectance_Band1','Nadir_Reflectance_Band2',
        'Nadir_Reflectance_Band3','Nadir_Reflectance_Band4',
        'Nadir_Reflectance_Band5','Nadir_Reflectance_Band6',
        'Nadir_Reflectance_Band7']);
  var orig = image;
  var region = Extent;
  var scale = 500;
  var bandNames = ['Nadir_Reflectance_Band1','Nadir_Reflectance_Band2',
  'Nadir_Reflectance_Band3','Nadir_Reflectance_Band4','Nadir_Reflectance_Band5',
  'Nadir_Reflectance_Band6','Nadir_Reflectance_Band7'];
  var meanDict = image.reduceRegion({
    reducer: ee.Reducer.mean(),
    geometry: region,
    scale: scale,
    maxPixels: 1e9
  });
  var means = ee.Image.constant(meanDict.values(bandNames));
  var centered = image.subtract(means);
  var getNewBandNames = function(prefix) {
  var seq = ee.List.sequence(1, 7);
  return seq.map(function(b) {
    return ee.String(prefix).cat(ee.Number(b).int());
    });
  };
  // PCA function
  var getPrincipalComponents = function(centered, scale, region) {
    var arrays = centered.toArray();
    var covar = arrays.reduceRegion({
      reducer: ee.Reducer.centeredCovariance(),
      geometry: region,
      scale: scale,
      maxPixels: 1e9
    });
    var covarArray = ee.Array(covar.get('array'));
    var eigens = covarArray.eigen();
    var eigenValues = eigens.slice(1, 0, 1);
    var eigenVectors = eigens.slice(1, 1);
    var arrayImage = arrays.toArray(1);
    var principalComponents = ee.Image(eigenVectors).matrixMultiply(arrayImage);
    var sdImage = ee.Image(eigenValues.sqrt())
    .arrayProject([0]).arrayFlatten([getNewBandNames('sd')]);
    return principalComponents.arrayProject([0])
    .arrayFlatten([getNewBandNames('pc')])
    .divide(sdImage);
    };
  var pcImage = getPrincipalComponents(centered, scale, region);
  return ee.Image(image.addBands(pcImage));
};

var modisColClip = modisCol
.map(function(image){return(image.clip(Extent))})
.map(filterqa)
.map(PCA)
.sort('system:time_start');
print("modis",modisColClip);

This code produces a error due to the masking of pixels. How can I change this? I have tried to swap around the order of functions to map, but I don't want to conduct a PCA using QA bands, and if I ignore the QA bands, the filter function fails because there is no QA bands.
Any idea?


texts: How to add Legend Panel to Dojo / OpenLayers?

texts:
I have OpenLayers app on Dojo framework. I've been serving WMS layers with GeoServer and wanted to use GetLegendGraphic to build legend for all layers.

I've seen there are some ready widgets (like Legend Panel in Esri API or GeoExt) and was wondering is there anything I can easily use with Dojo or OpenLayers?


texts:
I'm trying to change the active layer in ArcIMS from the default top position to index 2 in the TOC.  The TOC has 3 grouped layers, but it is the first grouped layer that I want the active layer to be updated to position 2 not the default position 0.  I found this code link from esri showing how to update the active layer, however I can not seem to figure it out.  Here is the code from the aimsLayers.js file.

//ActiveLayerIndex=<number of the preferred active layer>;
    ActiveLayerIndex=1;

    if (ActiveLayerIndex>=layerCount) ActiveLayerIndex = 0;
    if (!LayerIsFeature[ActiveLayerIndex]) {
        var chk = 0;

        for (var i=layerCount-1;i>=0;i--) {
            if (LayerIsFeature[i]) chk = i;
        }
        ActiveLayerIndex = chk;
    }

I added the "ActiveLayerIndex=" and changed it to "ActiveLayerIndex=1;" the active layer moved one down to position 1 in the TOC.  When I change it to 2 it did not move to position 2 in the TOC and displays no active layer.  When I change the index in this line, if (ActiveLayerIndex>=layerCount) ActiveLayerIndex = 0;, it does not update the active layer either.  How can I get the active layer to be on the third layer (position 2) in the TOC?

UPDATE

I found another code snippet and was wondering by using it to point to the specific layer name in the TOC if it will re-assign the active layer?  Looks like it could be called within the TOC.htm

//global reference to mapframe
    var t = parent.MapFrame;

    function layerName(theFormRef) {
        var theLayerNumber = -1;
        var theLayer = theFormRef.value;

        for (var i=0;i<t.LayerID.length;i++) {
            if (t.LayerName[i] == theLayer){
                theLayerNumber = i
            }
        }
        if (theFormRef.checked) {
            t.LayerVisible[theLayerNumber] = 1;
            //activating layer, NB! raster layers needs to be excluded
            setActiveLayer(theLayerNumber);
        } else {
            t.LayerVisible[theLayerNumber] = 0;
        }
        t.sendMapXML();
    }

http://forums.esri.com/Thread.asp?c=64&f=786&t=278849


texts: How to get RGB fill color values from polygones?

texts:
I'm trying to get color values from layers inserted into canvas. I used this code but it seems that I'm doing something wrong:

for layer in QgsMapLayerRegistry.instance().mapLayers().values():
    sym = layer.rendererV2().symbols()[0]
    fill = QgsSimpleFillSymbolLayerV2().fillColor()
    p=fill.getRgb()
    print(p)

I always get (0,0,255,255) as a result.

Also this is how it looks Layer properties from one of the polygones in canvas.

 Is there any way or code to do this right?

I would also like to add that second fill color is added  using this code:

sym = layer.rendererV2().symbols()[0]
fill = QgsSimpleFillSymbolLayerV2()
fill.setFillColor(QColor.fromRgb(result[0], result[1], result[2]))
sym.appendSymbolLayer(fill)


texts: QGIS 2.18 selected features not working in processing algorithms

texts:
After reading up on processing algorithms I see there is not an option to have a "selected features" checkbox.  It states in the documentation:

Algorithms follow the following simple rules when using a vector layer.

If the layer has a selection, only selected features are used.
If there is no selection, all features are used.

However, when I run a process such as 'Random selection' from the Processing Toolbox all items (in my case all points), are processed, not just the features I have selected.  I have gone into the Processing - Options - General menu and have selected/deselected 'Use only selected features' but no success?

Am I not proceeding correctly or is/was there a bug in v2.18?


texts: GeoServer sub-pages and OWS service only for those logged-in

texts:
I am trying to edit source files of GeoServer to meet my project requirements. I want to make newly added OWS services (like "org.geoserver.hello" example) and some of the sub-pages (like "helloweb" example) accessible only for logged-in users.

Is there an example, source file or documentation page (apart from http://docs.geoserver.org/2.7.x/en/user/security/index.html and http://docs.geoserver.org/latest/en/developer/programming-guide/security/index.html) I could rely on?


texts:
I'm trying to use r.reclass fonction (grass7) on QGIS (2.18.15) for Mac but when I ran it, there was an error :

[Errno 13] Permission denied: u'tmp152292116024' See log for more details

I searched to solve it but I didn't find anything...
Somebody can help me ?


texts:
I am trying to move a vertex in QGIS 3.2. I can select nodes in QGIS 3.2 (which show up blue), but every time I try to move it with the node tool it just deselects the vertex.

How do I move a node in QGIS 3.2?


texts: Using st_distance to calculate shortest distance between a table with point geometry to a table with multipolygon geometry

texts:
I'm trying and failing to do a calculation of the shortest distance for each point in a point geometry table (15000 records) to the nearest polygon edge in another table (170000 records). The end result I want is a distance for each point record to the nearest polygon record. My tables looks like this:

CREATE TABLE "brondby_adresserV2"(
  kvhx character varying,
  coordx integer,
  coordy integer,
  pk serial NOT NULL,
  the_geom geometry(Point,25832),
  dist_coast double precision,
  dist_lake double precision,
  CONSTRAINT "brondby_adresserV2_pkey" PRIMARY KEY (pk)
)

and this:

CREATE TABLE soe
(
  gid serial NOT NULL,
  fot_id double precision,
  geom geometry(MultiPolygon,25832),
  CONSTRAINT soe_pkey PRIMARY KEY (gid)
)

I've tried using John Barca's answer from this post: PostGIS nearest point with ST_Distance where he uses the st_dwithin so that the query uses the geometry index.

The result I get from my query pussles me and something is definately wrong. When visualising the result in qgis, it looks like the calculation is not done to all the polygons.
This is the code that I've used:

update "brondby_adresserV2" set dist_lake = t.foo
from
(
select distinct on ("brondby_adresserV2".kvhx) "brondby_adresserV2".kvhx,
"soe".fot_id, st_distance("brondby_adresserV2".the_geom, "soe".geom) as foo
from "brondby_adresserV2", "soe"
where ST_DWithin("brondby_adresserV2".the_geom, "soe".geom, 5000)
order by "brondby_adresserV2".kvhx, "soe".fot_id,
st_distance("brondby_adresserV2".the_geom, "soe".geom)
) t
where "brondby_adresserV2".kvhx = t.kvhx

Can anyone guide me as to what I'm doing wrong?


texts: What does the default crs being EPSG:3857 in Leaflet mean?

texts:
In the documentation of leaflet, it says the following about the projection that is being used:

  L.CRS.EPSG3857  The most common CRS for online maps, used by almost
  all free and commercial tile providers. Uses Spherical Mercator
  projection. Set in by default in Map's crs option.

(emphasis mine)

From this I concluded that I should convert my maps to EPSG:3857 so that this can work in Leaflet. However, this gives very large numbers as coordinates. Which in fact does not work. Turns out that EPSG:3857 is not used in Leaflet. I learned this from here and here. Instead, EPSG:4326 is used.

My question is: what does this part in the documentation mean? I feel like it may need to be changed, but perhaps I don't understand what the map's crs options are...


texts: AttributeError: 'int' object has no attribute 'save' pops up in save folder

texts:
I am facing this error time and again when I am try to read my img rasters and running a conditional statement. I am new to Python. The code ran for one raster that is without for loop but not later. Here is the code:

import arcpy, os
from arcpy import env
from arcpy.sa import *
env.overwriteOutput = True
arcpy.env.workspace = r"C:\Users\libpub\Desktop\try"
arcpy.CheckOutExtension('Spatial')

rasterlist = arcpy.ListRasters()
for raster in rasterlist:
    out1 = Con(('Raster')>0.2, 1, 0)
    print out1
    out1.save("C:\Users\libpub\Desktop\try")

print "done!"


texts: is there any possibility to save style of layers and pop-ups designed in Portal for ArcGIS? If not, what are your best way to apply same style for layers and their pop-ups for new-published layer? Is there any clever shortcut?


texts: Assigning incrementing id for each feature group in attributes with QGIS expressions

texts:
I have a layer with a repeating field name, thus different groups of features (I also have a numerical group field). How can I add a field group_id, using QGIS expressions? What I want is kind of a uniqe identifier per group, counting from 1 to the last element contained in the group (see screenshot).
What I tried: Using array_length (array_agg( "group", "group")), I get the number of features per group (field no_in_group). With  generate_series (1, [features per group]), I get an array, counting from 1 to the last no. of each group.
However, I'm stuck how to assign each element of the array to one feature of the group. I suppose that array_foreach () will do the job, but I'm unsure how to use it in this case.
Screenshot: the field with red outline is what I want to get: a unique, incrementing id, starting from 1 for each group. Here, I introduced the values manually to make clear what I'm looking for:


texts: django "GeometryField" widget type to accept coordinates in floating point format

texts:
I'm trying to make a django form for user to specify the coordinates on a map and i did this by including small map widget where it's possible to choose geometry type as well as put it on map. But i'm looking for some widget where i can choose geometry type (like LineString, Point, Polygon) and specify values in separate field. Is there any existing solutions for this or i should add my own custom fields?

Now i have this form field:

geometry = forms.GeometryField(required=False)


texts: Segmentize very long lines in Openlayers

texts:
Is there a Segmentization method for lines in Openlayers?

My problem is when drawing a very long line in Openlayers it draws it as a straight line, although it should appear as a curve (because of the curvature of the earth)

See the attached figure (the figure is from SQLSever Management Studio), it's for the following line:

'LINESTRING(26.16867957625264 25.30479198505643, 126.3173281947216 25.30479198505643)'

In PostGIS there is ST_Segmentize that does exactly that, but I need to do that in the front-end.


texts: Extracting raster values to polygon: error "arguments imply differing number of rows:1549, 0"

texts:
I am trying to extract the mean of a raster data to polygons with the following code:

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

agri_suitability=raster("suit")
regions <- readShapeSpatial("Regions.shp") #SpatialPolygonsDataFrame

plot(agri_suitability)
plot(regions, add=T)

polygon_agri_suitability=extract(agri_suitability,regions,fun=mean,na.rm=T)
regions@data = data.frame(regions@data, polygon_agri_suitability)

write.csv(regions,"suitability_mean.csv")

After data.frame I get an error:

  arguments imply differing number of rows:1549, 0

My "Regions.shp" file in fact includes 1549 rows. But I don't know what this error means and how I can solve it.


texts: How can I get a confusion matrix using specific columns in two different rasters in ArcGIS ? I have one raster of risk classes 1,2,3,4 and another raster of land cover A, B, C,D. I am looking for a matrix showing number of cell count for each land cover in each risk class.


texts: Can't run mapproxy server: SyntaxError: invalid syntax

texts:
I have successfully installed a MapProxy == 1.8.0. And I have set mapproxy-util create -t base-config gasps
But when I run server, using mapproxy-util serve-develop mapproxy.yaml

Console throws it:

l
[2019-09-04 11:36:08,452] mapproxy.config - INFO - reading: /root/mapproxy/gasps/mapproxy.yaml
Traceback (most recent call last):
  File "/root/mapproxy/bin/mapproxy-util", line 10, in <module>
    sys.exit(main())
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/script/util.py", line 365, in main
    commands[command]['func'](args)
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/script/util.py", line 84, in serve_develop_command
    app = make_wsgi_app(mapproxy_conf, debug=options.debug)
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/wsgiapp.py", line 106, in make_wsgi_app
    services = conf.configured_services()
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/config/loader.py", line 209, in configured_services
    return self.services.services()
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/config/loader.py", line 1479, in services
    new_services = creator(service_conf or {})
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/config/loader.py", line 1627, in demo_service
    layers[layer_name] = layer_conf.wms_layer()
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/util/py.py", line 78, in wrapper
    cache[args] = func(self, *args)
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/config/loader.py", line 1326, in wms_layer
    from mapproxy.service.wms import WMSLayer
  File "/root/mapproxy/lib/python3.7/site-packages/mapproxy/service/wms.py", line 36
    from mapproxy.util import async
                                  ^
SyntaxError: invalid syntax

asyncio is installed. How to solve a problem?


texts:
I received a shapefile for my project from an agency and was told that they always use UTM. However, the values have too many digits for being projected in UTM.

      min      max
 x 32507319 32624343
 y  5925974  6074881

I've tried reading the file into QGIS but the shapes always appear outside of my map's boundaries, even if I use a world map, no matter which projection I try. The area where they should appear is Schleswig-Holstein, which is the northernmost federal state of Germany and has the following min and max values in UTM

    min       max
x  426167.5  650147.9
y 5913377.7 6101334.8

As one can easily see, something's very wrong with the upper data. Does anyone know which projection I could use to fix this?


texts: Advanced python labeling in ArcMap - working with null values

texts:
I'm trying to label certain features using a specific field (osm_name_58_en).  However, when this field is null, I want to label the feature with (som_english_32_name).  If both are null then I want to label with the (amenity) field.

This is the code I'm trying but I'm not getting anything in return.  Maybe I'm completely out to left field? I'm a python and advanced labeling newbie.

def FindLabel ( [amenity], [osm_english_32_name], [osm_name_58_en] ):
    if not ([osm_name_58_en] is None):
        return [osm_name_58_en]
    elif ([osm_name_58_en] is None and not([osm_english_32_name] is None)):
        return [osm_english_32_name] + '\n' + [amenity]
    else:
        return [amenity];


texts:
I created a map with some marker which can be selected and have them name appears in an external field. But i even if i am able to select multiple marker on the map ( shift command ), only one name appears in my field.

I created a vector layer loading geojson content and an event with a "featureselected" on this layer.
But still nothing, here is my code.

var input = document.getElementById('input-airports');

var layerListeners = {
    featureselected:   function(e) {
        var value = '';
        value = e.feature.attributes.code  + ', ' + ' ';
        input.value = value;
    }
};

var vector = new OpenLayers.Layer.Vector("Airport Markers", {
    projection: "EPSG:4326",
    eventListeners: layerListeners,
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: new OpenLayers.Protocol.HTTP({
        url: "data.geojson",
        format: new OpenLayers.Format.GeoJSON()
    }),
    "select": new OpenLayers.Style({
       pointRadius:20
    })
});
map.addLayer(vector);


texts:
I would like to select this row (it is a polygon) from this fusion table (called CDM_areas). I am struggling with the code to select such a field...

I have so far (uncertain in bold):

var ndviclip = ndvichange.clip(CDM_areas.select('name' **= '40'**));


texts:
I am trying to perform a reprojection for data with datum Nad27 to WGS84 datum. According to AndreJ in some previous questions in Stackexchange (CRS reprojection for instance) QGIS uses internally by default the ntv2 transformation from NAD27 to WGS84, that should be very exact. I tested with the custom crs calculation tool. First a lat lon point at 99w , 19n with parameters for EPSG32614 looks like this:

The same point (which is near Mexico City) with EPSG 26714 parameters looks like this:

There is a difference of 126 meters on the north coordinate but no difference in the east coordinate which is odd. I tested a custom crs I created with a towgs84 argument as follows:

Now I have a difference for both the east and north coordinate that should be expected for the datum change. Adding the same towgs84 to EPSG 26714 makes no difference with figure 2

Making the reprojection of the same coordinates in GPS Babel  gives me that the point in UTM NAD27  is 500032mE, 2,100,627mN. There is almost no difference with my custom crs. So I wonder how to make the correctly the reprojection of the data.


texts: Determining which polygons line crosses and length of it within each polygon using ArcGIS Desktop?

texts:
I have two layers, one with polylines representing cables in the ground and one layer of the land cover. The cable data consists of around 200 different objects, all with an ID and a total length. I would like to get to know which types of land cover each cable cuts through and I need it to be connected to the line's ID.

I've tried to do an intersect with the line layer and the land cover layer but the problem then is that when one line cuts through more than one type of land cover I get several different polylines for that feature. I would prefer if I could get the different land cover types as fields in the attribute table so that I can see how much each cable crosses each land type


texts: Python code for transforming lat long into FIPS codes?

texts: I am trying to transform CSV files of lat long into FIPS codes. Does any one have a Python code that is able to do this? I have that worked then stopped working giving me an 'http 400 bad request' output.


texts: Multiple polygon to shapefile in XYZ format


texts:
Are you allowed to develop closed source, commercial plugins for QGIS?
License is GNU GPL, which is normally very restrictive.

Update: Assume it is a python plugin (possibly open-sourced), that call's a closed-source dll, which doesn't link to anything from QGIS.


texts: Display points when there's great distance between them (GWT-Openlayers)

texts:
The case is the following:
I have a layer and there are two points on it. The first is in Australia, the second is in the USA. The continent or the exact position of the points doesn't count. The essential part is the great distance between the points.
When the application starts, the first point appears (zoomlevel is 18). The second point isn't displayed because it is far away from here and the zoomlevel is high. Then i call the panTo function with the location of the second point and the zoomTo function with the same (18) zoomlevel. The map jumps to the right location but the second point doesn't appear. The point appears only if i zoom in/out or resize the browser window.

I assume Openlayers doesn't load the points which aren't on the screen actually. And those points are loaded only when i jump to their location (panTo) and then change the zoomlevel (or resize the browser window).
Has anybody experienced such behaviour? Or has anybody got a solution to this problem?

Thanks in advance!


texts:
(I am using QGIS 2.6 Brighton)

I have 2 shps, one points and another one polygons. I need to make them one.

Is there a way to join them or convert the points to polygons? Vonoroi doesnt work for me this time as the points are individual plant locations. I thought that might be possible to create regular polygons (1mx1m) instead of each point. Ore anyother size, really....

I've tried "points2one", I thought that every point could be a different little polygon... but didn't work or I didn't know how to do it (if this is possible....)

Any ideas?


texts: Why are my icosahedron triangle subdivisions not equilateral (PostGIS)?

texts:
I am creating a geodesic polyhedron using PostGIS. I am dividing the base icosahedron points (see Appendix A), into the first level of Class I subdivisions.
The setup
For simplicity let's take just the first triangle formed by these points:
CREATE TABLE pix (
    id bigserial
  , name text
  , geog geography(POINT,4326) -- use WGS84
);

insert into pix values (DEFAULT, 'China', st_point(122.3, 39.1));
insert into pix values (DEFAULT, 'Norway', st_point(10.53619898, 64.7));
insert into pix values (DEFAULT, 'Arabian sea', st_point(58.15770555, 10.44734504));

Which creates the base icosahedron triangle

The numbers of the triangle sides correspond to the row numbers shown in the table generated by this query
select p1.name
  , p2.name
  , st_distance(p1.geog, p2.geog) as dist
from pix p1
  cross join pix p2
where p1.id < p2.id

The points are roughly equidistant. There are small variances because the seed data above seems to have been generated assuming a sphere, not spheroid, which st_distance is returning - use_spheroid false in st_distance will yield much closer distances, but note that using false in all the presented queries does not change the problem described below.
First level of Class I subdivision
insert into pix (name, geog)
select p1.name || '-' || p2.name
  , st_project(p1.geog, st_distance(p1.geog, p2.geog) / 2, st_azimuth(p1.geog, p2.geog)) as geog
from pix p1
  cross join pix p2
where p1.id < p2.id

The above generates the mid-points of each of the pairs of points on the original triangle

And these have the following distances, generated by this query
select p1.name
  , p2.name
  , st_distance(p1.geog, p2.geog) as dist
from pix p1
  cross join pix p2
where p1.id < p2.id
  and p2.id > 3
order by dist

The problem
The distances of the segments 7, 8 and 9 (around 4000kms) are much longer than the distances of segments 1-6 (around 3500kms).
Why is this algorithm not creating something more close to equliateral triangles? Is my use of st_azimuth and st_project wrong? Or did I miss some more basic maths about projections of points onto spheres?


texts:
ArcScene includes the capacity to animate 3D text added to the scene, through the use of Graphics Layers. I'm able to use this to make text fade in and out as visibility of layers is changed, however I'd like to achieve something similar in ArcMap.

I've got a clunky work around which involves creating new point layers with transparent symbology, labeling them to suit the different layers, and fading those in with the relevant layer. I can't apply the labels to the layers I'm displaying as they are all raster. Is there a more straight forward way to achieve this?

I could do this in ArcScene with a view point set perpendicular to the plane of the data, but I'd rather use ArcMap because it is more stable with the amount of raster data I'm using.


texts: How to select a fixed dimension ROI in ENVI?

texts: ENVI provides the freedom to select ROIs by drawing polygons etc. I want to clip a fixed size ROI, say 512 x 512 in size. Is there a way to do that?


texts: Copy/Paste Qgis Plugin or derive new Plugin from existing one

texts: I created a Python Plugin for QGIS and I want to experiment with various versions of this plugin. Is there an easy way to copy/paste a plugin or to create a plugin from another plugin? ( Otherwise one should have to repeat manually the steps of developing the plugin, or copy/paste the plugin folder and rename the file names and all include statements)


texts: ArcMap Calculate Average Nearest Neighbor - failed

texts:
I am trying to determine the mean and standard deviation from a subset of point features. Once I determine the standard deviation I want to exclude any points that are more than twice the standard deviation from the mean distance.

To do this I tried to run Average Nearest Neighbor for a selected set of records in a WGS_1984 point file using ArcMap 10.2.2 basic. I've excluded ungeocoded records from the analysis and checked to see if my study area was greater than 30 degrees. Keep getting this error:

Traceback (most recent call last):
  File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\NearestNeighbor.py", line 348, in <module>
    setupNearestNeighbor()
  File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\NearestNeighbor.py", line 55, in setupNearestNeighbor
    nn = NearestNeighbor(ssdo, concept = concept, studyArea = studyArea)
  File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\NearestNeighbor.py", line 106, in __init__
    self.setStudyArea()
  File "c:\program files (x86)\arcgis\desktop10.2\ArcToolbox\Scripts\NearestNeighbor.py", line 124, in setStudyArea
    spatialRef = ssdo.spatialRef)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Scripts\SSUtilities.py", line 442, in innerFunction
    returnValue = functionName(*args, **kw)
  File "C:\Program Files (x86)\ArcGIS\Desktop10.2\ArcToolbox\Scripts\SSUtilities.py", line 1910, in minBoundGeomPoints
    "", "MBG_FIELDS")
  File "c:\program files (x86)\arcgis\desktop10.2\arcpy\arcpy\management.py", line 2724, in MinimumBoundingGeometry
    raise e
ExecuteError: ERROR 000210: Cannot create output :\ArcGIS\Default.gdb\regularBound_FC0
Failed to execute (MinimumBoundingGeometry).

Thought maybe it was the length of the output name or a lock on the default gdb but get same results even after changing default gdb. Can't manually change the output file as you do not have this option in the tool parameters.

Would like to know if anybody else has run into this problem or knows of a better process for calculating distance standard deviations and means.


texts: How to find the proper coordinate reference system with two sets of data using different measurements or too convert one system into another

texts: I have two different entries of data in two different spread sheets. One of the spread sheets uses a coordinate reference system that I have no idea what the units of measurement are. the properties of the shape file say the crs is +proj=longlat +datum=WGS84 +no_defs, however the data points show up off in the middle of nowhere. The other set of data is properly aligned to my map and using the same crs. Both sets of data were collected from a hand held gps which i no longer have access to. Is there some way for me to convert my data all at once to realign properly with my map? Or do I have to try using a different crs? Finally if I do use a new CRS will that not allow me to merge the layers or do types of statistical analysis?


texts: I have created specific fields for 'Creator', 'Create Date', 'Editor' and 'Date Edited' and have enabled Editor Tracking on these fields. Whenever I add an attribute, my username and the date are populated in all four fields. Why is this happening? It should only be populating my 'Creator' and 'Create Date' upon feature creation, and only populate the 'Editor' and 'Date Edited' once that feature gets edited, correct?


texts: Illinois Legislative Zip to District Match Data Set

texts:
I have the Zip +4 for the about 650 areas and I need to find Illinois Senator and Rep code,
address, name. How can I solve this issue. Did anyone solve this issue using any vendor data? We do not want to look up the names of the senators from the website.


texts:
I been using mission planner for some aerial survey study (as my ground control station and autonomous mission planner). However, I haven't been able to find all the information to properly cite this software.

Mission Planer; Computer Software; ArduPilot Development Team: "City, Country", 2012.


texts: Seamlessly combine map layers (fill the gaps of one layer using the second one)

texts:
I am creating air quality maps using two approaches.

One of them gives slightly better results when validated against measurements but there are significant spatial gaps present in the output.

The other approach gives a spatially continuous layer.

My goal is to combine the two layers by taking as much information as possible from the first layer and fill in the missing data with information from the second layer. I would like to avoid simply averaging the layers as this 1) does not retain the information from the more "valuable" layer, 2) leads to sharp breaks between the layers at some of the edges (as shown in the last picture where unnatural squares appear).

Is there any more sophisticated method to combine data in this way?

For example by assigning weights to the border data to combine it gradually.

Ultimately I would like to accomplish this in R.


texts:
Is there a way in ArcGIS 10.1 to have a text field with a defined width to perform word wrapping (ideally with hyphenation) automatically, or do i have to implement it myself using e.g. pyhyphen?

Use case are long (>200 character) attribute values that need to be displayed. At the moment i do hard linebreaks after every 50'th character, but that is not pretty.


texts: Deactivating measurement tool in OpenLayers 5?

texts:
I have created a web portal in which measurement tool is available where user select measure distance and area.

Now my problem is how to activate and deactivate it when user click on button.
here is my html

 <form class="measuring-types">

            <input type="radio" name="measurement" value="Polygon" checked> Area

            <input type="radio" name="measurement" value="LineString"> Distance

    </form>

here is my javacsript code

var measurementRadios = $('[type=radio]');
var resultElement = $('#js-result');
var measuringTool;

var enableMeasuringTool = function() {
  map.removeInteraction(measuringTool);

  var geometryType = measurementRadios.filter(':checked').val();
  var html = geometryType === 'Polygon' ? '2' : '';

  measuringTool = new ol.interaction.Draw({
    type: geometryType,
    source: vectorLayer.getSource()
  });

  measuringTool.on('drawstart', function(event) {
    vectorLayer.getSource().clear();

    event.feature.on('change', function(event) {
      var measurement = geometryType === 'Polygon' ? event.target.getGeometry().getArea() : event.target.getGeometry().getLength();

      var measurementFormatted = measurement > 100 ? (measurement / 1000).toFixed(2) + 'km' : measurement.toFixed(2) + 'm';

      resultElement.html(measurementFormatted + html);
    });
  });

  map.addInteraction(measuringTool);
};

measurementRadios.on('change', enableMeasuringTool);

enableMeasuringTool();

How to deactivate measurement on button click.


texts:
I need to create a table based on a Python list and paste it in layout view in QGIS.

I also found a QgsLayoutTable class in PyQGIS but there's no example of creating a simple table.

Does anyone know how to perform that?


texts: Mapping application basic for ruby on rails

texts:
I am new to build map application using ruby on rails, with rails 4.1.1.

I want to create map using map tiles. i have to use PostgreSQL database with postgis.

Can you please explain how this flow will work ? Form where i have to start for map creation and display.

Thanks,

Ghanshyam


texts:
I am trying to cut a window (in pixels) from a geotiff file and output in a different size. gdal_translate allows that but gives no option to set resampling, while gdalwarp allows resampling while it requires the src window to be in georeferenced coords.

Is there away to get the best of both. i.e. (src in pixels and using resampling)?

Thanks.


texts: I want to add an alert message to my python add-in(which will be like, when I click a button it will summon a warning message ). But I tried error message functions like : AddError, AddWarning and they just dont do anything.


texts:
I have osm.xml file where projection is defined as follows:
Map background-color="grey" srs="&srs900913;" minimum-version="2.0.0">

n settings.xml file projection is also set there as:

<!ENTITY osm2pgsql_projection "&srs900913;">

I want to know what is the difference between these 2 projection. I have changed osm2pgsql projection to wgs-84 and it runs f9 for some files. i am confused in this thing

EDIT

When i use srs4326 projection i am getting error on console which is as follow: Unknown projection string, using web mercator as never the less. +init=epsg:4326

I have long/lat coordinates i want to use srs4326 projection


texts: How to Divide geo spatial Raster images into tiles and combine back into a complete image

texts:
I am reading a 2048 x 2048 image in the form of (204 x204) tiles and bilinear interpolating  them to get (819x819) tiles. Then combining several 819x819 tiles to form a 8190x8190 image.
Now as you can see in the attached image, after interpolation the image has tiled edges, Is there anything wrong with interpolation or with the way the tile is being written.

for (j = 0; j < nYsize_R; j = j + nYsize_R_1)
    {

    if (j == nYsize_R - (nYsize_R_1 + nYsize_R_2))
    {

        nYsize_R_1 = nYsize_R_1 + nYsize_R_2;
    }

    i2 = 0; ii = nXsize_P_1; nXsize_R_1 = temp_nXsize_R_1;
    for (i = 0; i < nXsize_R; i = i + nXsize_R_1)
    {
        if (i == nXsize_R - (nXsize_R_1 + nXsize_R_2))
        {

            nXsize_R_1 = nXsize_R_1 + nXsize_R_2;

        }
        if (i2 == nXsize_P - (nXsize_P_1 + nXsize_P_2))
        {
            ii = ii + nXsize_P_2;

        }

        poBand_MS[0]->RasterIO(GF_Read, i, j, nXsize_R_1, nYsize_R_1, poBandBlock[0], nXsize_R_1, nYsize_R_1, GDT_Float32, 0, 0)

        poBandBlockB2 = resizeBilinear(poBandBlock_padded[0], nXsize_R_1 + 1, nYsize_R_1 + 1, nXsize_P_1, nYsize_P_1);

        poBand_MS_interpolated[0]->RasterIO(GF_Write, i2, j2, ii, jj, poBandBlockB2, ii, jj, GDT_Float32, 0, 0);

        i2 = i2 + ii;

    }
    j2 = j2 + jj;
    if (j2 == nYsize_P - (nYsize_P_1 + nYsize_P_2))
    {
        jj = jj + nYsize_P_2;
    }

}

float* resizeBilinear(float * pixels, int w, int h, int w2, int h2) {

double a, b, c, d;
int x, y;
unsigned int    index;
float x_ratio = ((float)h-1) / h2;
float y_ratio = ((float)w-1) / w2;
float x_diff, y_diff, blue, red, green;
int offset = 0;
int  xx = 0;

for (int i = 0; i<w2; i++) {
    for (int j = 0; j<h2; j++) {
        x = (float)(x_ratio * j);
        y = (float)(y_ratio * i);

        x_diff = (x_ratio * j) - x;
        y_diff = (y_ratio * i) - y;
        index = (y*(h) + x);
        a = pixels[index] ;
        b = (unsigned int)pixels[index + 1];
        c = (unsigned int)pixels[index + (h)];
        d = (unsigned int)pixels[index +( h) + 1] ;

        Buffer2[xx] = (unsigned int)(a *(1 - x_diff)*(1 - y_diff) + b * (x_diff)*(1 - y_diff) + c * (y_diff)*(1 - x_diff) + d * (y_diff)*(x_diff));
        //cout << Buffer2[xx] << endl;
        xx = xx + 1;

    }
}
return Buffer2;

}


texts: Using data captured through a GP task window in a Flexviewer Print Widget

texts:
I have a Flex Viewer that allows users to undertake a GP raster processing task by selecting a series of datasets (rasters with predefined scientific scores) through the GUI and then asks the user to assign an inputted 'weight' (1-3) score for each variable before the analysis is undertaken.

The 'weights' are already captured as parameters in my model and used in the raster calculation process ("%Dataset1%" * float(%Theme1 Weight%)).

I have the mapped output displaying in the print window and would like to also show the 'weights' that were inputted by the user (Theme1 Weight).

Any suggestions on how this could be done?

The image below is the GP input window for users. What I want to do is display the Weight scores input by the user (2 for Biodiversity and 1 for Air etc) in the print window.


texts: Editing KML custom data elements using Google Earth?

texts:
I've been reading up a bit on KML extended tags and custom data elements.  https://developers.google.com/kml/documentation/extendeddata

I have a person that likes using Google Earth but he is also in charge of updating attributes (and some features) but he refuses to use a GIS like QGIS or ArcGIS.

I can get a kml file formatted with all the attributed he needs for each feature, using the extended data tags, but I want him to be able to change them and then give me the file back so it can be integrated into our GIS.

Does anyone know of an interactive editor that allows editing of "attributes" formatted as extended data elements using Google Earth.


texts: How to Split a vector file into several different files using a Grid in QGIS?

texts:
I have a shapefile containing contour for a large area, and I want to split it into different shapefiles each covering a 1KM X 1KM area.
I have already created the required Grid shapefile in QGIS.  How do I now split the contours using the created Grid?

I am searching for something equivalent to the Split Tool in ArcGIS.


texts:
I am new to ArcObjects, I need to identify the angle of the line in a particular Layer "MV Feeder"

I am trying the below code, which gives error "object reference is not set to and instance of the objects":

 public void ShowLineAngle()
    {
        IMxDocument mxdoc = ArcMap.Application.Document as IMxDocument;
        IMap map = mxdoc.FocusMap;
        if (map.Layer[0] == null)
        { return; }

        for (int i = 0; i < map.LayerCount; i++)
            {
                if (map.Layer[i].Name == "MV Feeder")
                {
                    ILayer HTOverheadLine = map.Layer[i];

                    IFeatureLayer fl1 = new FeatureLayerClass();
                    fl1 = HTOverheadLine as IFeatureLayer;
                    ILine line = new LineClass();
                    IQueryFilter qF = new QueryFilterClass();
                    qF.WhereClause = "";
                    IFeatureCursor featureCursor = fl1.Search(qF, false);
                    IFeature feature = null;
                    while ((feature = featureCursor.NextFeature()) != null)
                    {
                        line = feature.Shape as ILine;
                        double angle = line.Angle;
                        System.Windows.Forms.MessageBox.Show("Angle = " + angle);
                    }
                }
            }


texts: Combing multiple LAS files with different vertical and horizontal accuracy?

texts:
The data set in question was made using .las and .zlas files from respectable sources.

Data set 1 was collected in the Fall of 2014. The RAW flight line swath data was processed to create classified LAS 1.2 files formatted to 742 individual 1500m x 1500m tiles. The vertical accuracy is 18.3cm at the 95 percent confidence level. The horizontal accuracy is 1.04m at the 95 percent confidence level. The comprehensive metadata files will be made available.

Data set 2 was collected January 25, 2012 through April 11, 2012. The tiles were formatted to 5,000m x 5,000m in ASPRS ver 1.2 format. The vertical accuracy is 0.36m at the 95 percent confidence level. The horizontal accuracy is 1.04m at the 95 percent confidence level. The comprehensive metadata files will be made available.

GIS workflow:

The data was processed using Arcmap 10.5.

The two datasets had different sized tiles that overlapped in several areas. The larger tiles below are from the NRCS 2012 dataset and the smaller tiles are the SDGS 2014 dataset.

The Tile LAS tool in the 3D Analyst toolbox was used to create a new LAS dataset with uniform tiles. Rearrange_Points option was checked under LAS options. This orders LAS points according to their spatial clustering. The target file size was set at 250MB.

The result was a new LAS dataset with uniform tiles 250 MB in size. This took about 7 hours of processing time.

The filter settings underneath the layer properties menu were changed to show ground points only.

After the LASD layer was changed to show ground points only, the LAS Dataset to Raster tool from the conversion toolbox was used create a raster dataset. This required a new file geodatabase to be created for storing the new output files. The interpolation type was set to triangulation and the method was set to natural neighbor. The grid cell size used was 10 x 10m. This took about 4 hours of processing time.

The resulting raster dataset is below and is now a bare earth DEM. The dataset came out as expected however there are some areas on the eastern part of the river that appear as white specs. I think these are most likely areas where there were no ground points.

Can this be fixed by using the Binning interpolation type?

My main question was to ask the forum if there is a better way to combine LAS files like this. I'd like to make another data set of the same DEM but use .las files from a third data set with higher vertical accuracy. I have reserves about this because I am not sure what the Tile LAS tool does to the data points. Like I don't know how using this tool affects the accuracy of the final DEM and if there is a way to keep track or recalculate vertical and horizontal accuracies for the final DEM. I have limited processing power and the final DEM will be used for stream geomorphology analysis.


texts: ST_ColorMap wrongly loses 0.0 value when INTERPOLATE


texts:
I have a raster layer that has cell size(x, y) 0.008928, 0.008928.
I download it from the website as 1km resolution data under WGS_84. I can not find the unit for these numbers in the metadata. Is there a proper way to calculate the cell size in meters? I need to resample them to 4km resolution data.


texts: Get minimum and maximum Latitude and Longitude of New York?

texts:
How can I find the minimum and maximum Latitude and Longitude of New York?

I need the coordinates of New York limits, in order to crop the area of the map outside of New York.


texts:
I am working on a city transport map, and for fun (and practice) I try to make an area of influence of transport station. But, whatever I try, something gets wrong. So two questions.

1 - How can I fix marker size with distance? I have metro, tram and bus stations. Each one has a different area size (1km, 750 ans 500m), that I have included in the datasheet already.

2- How can I fix marker size whatever the zoom is, and automatically speaking (that means I don't have to write in CSS code for each zoom level)?


texts: Raster image preview showing empty layer in geoserver

texts:
I have added a raster layer (GeoTIFF) on GeoServer, it is visible up to a certain zoom, but after that, it shows a blank screen, even though on click, it is returning values for RED, GREEN, BLUE, ALPHA band.

blank image

working image
Logs are showing no errors
2021-03-24 21:48:30,934 INFO [geoserver.wms] -
Request: getServiceInfo
2021-03-24 21:48:31,543 INFO [geoserver.wms] -
Request: getMap
    Angle = 0.0
    BaseUrl = http://localhost:8080/geoserver/
    Bbox = SRSEnvelope[344457.52820239705 : 344719.3898535189, 3117493.9743865044 : 3117952.083152021]
    BgColor = java.awt.Color[r=255,g=255,b=255]
    Buffer = 0
    Clip = null
    CQLFilter = null
    Crs = PROJCS["WGS 84 / UTM zone 43N",
  GEOGCS["WGS 84",
    DATUM["World Geodetic System 1984",
      SPHEROID["WGS 84", 6378137.0, 298.257223563, AUTHORITY["EPSG","7030"]],
      AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich", 0.0, AUTHORITY["EPSG","8901"]],
    UNIT["degree", 0.017453292519943295],
    AXIS["Geodetic longitude", EAST],
    AXIS["Geodetic latitude", NORTH],
    AUTHORITY["EPSG","4326"]],
  PROJECTION["Transverse_Mercator", AUTHORITY["EPSG","9807"]],
  PARAMETER["central_meridian", 75.0],
  PARAMETER["latitude_of_origin", 0.0],
  PARAMETER["scale_factor", 0.9996],
  PARAMETER["false_easting", 500000.0],
  PARAMETER["false_northing", 0.0],
  UNIT["m", 1.0],
  AXIS["Easting", EAST],
  AXIS["Northing", NORTH],
  AUTHORITY["EPSG","32643"]]
    Elevation = []
    Env = {}
    Exceptions = application/vnd.ogc.se_inimage
    FeatureId = null
    FeatureVersion = null
    Filter = null
    Format = image/jpeg
    FormatOptions = {}
    Get = true
    Height = 769
    Interpolations = []
    Layers = [org.geoserver.wms.MapLayerInfo@bc01a9e8]
    MaxFeatures = null
    Palette = null
    RawKvp = {EXCEPTIONS=application/vnd.ogc.se_inimage, FORMAT=image/jpeg, REQUEST=GetMap, SRS=EPSG:32643, BBOX=344457.52820239705,3117493.9743865044,344719.3898535189,3117952.083152021, VERSION=1.1.1, STYLES=, WIDTH=440, SERVICE=WMS, HEIGHT=769, TRANSPARENT=true, LAYERS=tiger:testgarud}
    RemoteOwsType = null
    RemoteOwsURL = null
    Request = GetMap
    RequestCharset = UTF-8
    ScaleMethod = null
    Sld = null
    SldBody = null
    SldVersion = null
    SortBy = null
    SortByArrays = null
    SRS = EPSG:32643
    StartIndex = null
    StyleBody = null
    StyleFormat = sld
    Styles = [StyleImpl[ name=raster]]
    StyleUrl = null
    StyleVersion = null
    Tiled = false
    TilesOrigin = null
    Time = []
    Transparent = true
    ValidateSchema = false
    Version = 1.1.1
    ViewParams = null
    Width = 440


texts: I have a new table with statistics data that has been calculated with the tool: Zonal Statistics as a Table; in order to know the min, max, mean, range, std, etc, from a thickness value that was previously calculated. My new table with this statistics data, has a field called ObjectID which is also in another old table, I want to relate this new statistic data with this old table in order to have all the information of each objectID (within both tables) in just one table. I am trying to use the Join tool, in order to do this, but maybe because there are almost 3000 ObjectID in these input table the computer keeps getting stack, I have also tried with the add join tool but it also doesn't work, what can I do?


texts: How to extract the speed from a .gpx file?

texts:
Sometimes the .gpx file that you obtain from your GPS-Device does not deliver the speed. The file looks like this:

[...]
<trkpt lat="59.4179482758" lon="3.6870684847">
    <ele>518.38</ele>
    <time>2016-07-17T08:27:18Z</time>
</trkpt>
[...]

The speed is missing.
How is it possible to obtain the speed?


texts: Auto-style shp layers with pre-defined schema in QGIS?

texts:
I have shapefiles with a pre-determined format, but different data.  Is there a way in QGIS to set the style automatically based on pre-determined fields, before any files are open?

I'm thinking something along the lines of a stylesheet which references the fields, possibly based on the file name.


texts: GeoJSON convert Polygons and MultiPolygons to LineStrings and MultiLineStrings

texts:
Is there a gdal or python based tool for converting Polygons and MultiPolygons inside of GeoJSON to LineStrings and MultiLineStrings? I'm using rasterio / rio rasterize or gdal_rasterize to rasterize GeoTIFF w/ GeoJSON.

I need to read in GeoJSON files and be able to convert any (Multi)Polygons into LineStrings and MultiLineStrings so I can rasterize just the boundary or outline of the (Multi)Polygons onto GeoTIFFs (this includes exterior and interior polygons) and not have it fill in the entire shape.

Any ideas?


texts:
I am new for GIS programming. I try to change geotiff coordinate system from epsg 32635 to epsg 4326. I used the following gdal command:

gdalwarp -rc -s_srs epsg:32635 -t_srs epsg:4326 -dstnodata 0 -of GTiff -co tiled=yes UU5L.tif UU5L_rej.tif

The source image is on the left and the result is on the right.

Obviously something went wrong. When I check metadata on result image, I see that corner coordinates in second cell are not the same as source image second cell.

What do these corner coordinates mean? I know first values are image place on Earth, but second?

Here is a little example:

Corner Coordinates:

Upper Left  (  500000.000, 7530000.000) ( 27d 0' 0.00"E, 67d53' 4.61"N)
Lower Left  (  500000.000, 7434000.000) ( 27d 0' 0.00"E, 67d 1'24.61"N)
Upper Right (  596000.000, 7530000.000) ( 29d17' 2.00"E, 67d52' 7.37"N)
Lower Right (  596000.000, 7434000.000) ( 29d12'10.21"E, 67d 0'29.74"N)
Center      (  548000.000, 7482000.000) ( 28d 7'18.05"E, 67d27' 0.64"N)

Now i'll understand more, today i face a new challenge.
I have two different geotiff files both has wgs84 crs (epsg:4326).

First tiff will show rectangular area, but second tiff looks like baseball field i mean top of map has curved line and bottom has short curved line, both side lines are also curved on the center of area.

So do i have any change to handle second tiff so it would looks like first one?

I have been using gdal scripts for handling, can i do that for gdal?


texts: WMS line layer is not visible in Openlayers

texts:
I have few shape files (Polygon/Lines/Points) imported into Postgres. I have published the Layers in Geoserver.

The Polygon layer is displayed properly in Openlayers However when I consume Line Layer in openlayers it is not visible neither is the point layer.

Also upon layer preview for line is not as easy as it was for polygon.Line layer preview is as below (it is 1 line , But 3 lines are displayed) :

I am adding Line layer in openlayers as below :

var layers = [
    new ol.layer.Tile({
        source: new ol.source.OSM()
    }),
    new ol.layer.Tile({
        //Polygon layers : Works Fine
        source: new ol.source.TileWMS({
                url: 'http://xxxxxxxxx:8080/geoserver/localhost/wms?',
                params: { 'LAYERS': 'localhost:village_maps', 'TILED': true },
                serverType: 'geoserver',
                // Countries have transparency, so do not fade tiles:
                transition: 0,
                transparent: true
                //format: new ol.format.GeoJSON()
            });
    }),
    new ol.layer.Tile({
        // Line Layer : Overlaps above polygon layer and Not visible
        // After hiding the above polygon layer also line layer is not visible.
        source: new ol.source.TileWMS({
            url: 'http://xxxxxxxxx:8080/geoserver/localhost/wms?',
            params: { 'LAYERS': 'localhost:proposed_line', 'TILED': true },
            serverType: 'geoserver',
            // Countries have transparency, so do not fade tiles:
            transition: 0,
            transparent: true
            //format: new ol.format.GeoJSON()
        });
    })
  ]

Please suggest.


texts: Checking if layer is visible in legend using ArcPy?

texts:
I want to be able to check if a layer is visible in the legend and if so then perform another operation - in python.  I am using ArcGIS 10.2 for Desktop.

I cannot find anywhere how to check if a layer is visible or not in the legend.


texts: NumPy is missing, but fixing the environment variable removed my printer settings?

texts:
I recently installed regular Python on my workstation in addition to ArcGIS 10.2. Based on various things I found online, it seemed like the consensus was to install the full Python, because there are some libraries that ArcGIS's version of Python can't access. I am a Python noob so am trying to learn as I go. I've surely messed something up, but I don't know how to fix it!

I have been getting variations of this error in ArcMap:

>>> import numpy
Runtime error
Traceback (most recent call last):
File "<string>", line 1, in <module>
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\__init__.py", line 24, in <module>
from arcpy.toolbox import *
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\toolbox.py", line 342, in <module>
from management import Graph, GraphTemplate
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\management.py", line 22, in <module>
import _management
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\_management.py", line 14, in <module>
import _graph
File "c:\program files (x86)\arcgis\desktop10.1\arcpy\arcpy\_graph.py", line 27, in <module>
import numpy
ImportError: No module named numpy`

I followed the directions at this link to set the path variable for Python to the specified library folder in Windows. It worked and I was able to use the tool that I needed.

When I opened the .mxd that I'm currently working on, which is an Arch E size, the layout board was letter-sized with the data frame being the Arch E size. When I checked under the print settings, the only page size shown for the plotter was letter size.

When I removed the environment variable, my map document was "fixed" and my paper sizes were back.

I don't have any other Python folders in the location that the article mentions, C:\Python27\ArcGIS10.2\Lib\site-packages. I'd like to be able to use the full install of Python, but if I have to, I can uninstall it. I didn't have this NumPy error before installing Python, so my guess is that I overwrote something I shouldn't have. How should I go about trying to fix this?


texts: Seeking Colombia municipality geometries

texts:
I'm looking for shapefiles (or some similar format) of Colombia's municipality boundaries.

Wikipedia has a map that seems to be correct and which matches another dataset that I will eventually be merging in based on municipality names and DANE codes. I've tried contacting the file creator about the underlying source data, but without success.

I've tried GADM, but some of the municipalities are missing/incorrect, e.g. all the way in the South Amazons should have 11 municipalities, not 8.

Here's OpenStreeMap with data export via Mapzen. Many municipalities are simply just missing, although the borders for those that are there appear to be correct and more detailed than GADM.

I guess I could try to combine departments that OpenStreeMap had all municipalities for with the GADM data, but it looks like that would still only partially get me there. And in general it just makes me a bit suspicious about using the GADM data at all.


texts: Are these two OGC WKT projections the same? And if so, why are their reference strings different?

texts:
So here is the .prj text of a federal shapefile dataset of all public airports in its native OGC WKT format:

GEOGCS["GCS_WGS_1984",DATUM["D_WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]]

Here is the OGC WKT text of a spatialreference.org match (http://spatialreference.org/ref/epsg/4326/ogcwkt/) which I think is effectively the same projection. Most of the time the two strings are an exact match though, so I'm not completely sure:

GEOGCS["WGS 84",DATUM["WGS_1984",SPHEROID["WGS 84",6378137,298.257223563,AUTHORITY["EPSG","7030"]],AUTHORITY["EPSG","6326"]],PRIMEM["Greenwich",0,AUTHORITY["EPSG","8901"]],UNIT["degree",0.01745329251994328,AUTHORITY["EPSG","9122"]],AUTHORITY["EPSG","4326"]]

As you can see, from a high level, these two look cannily similar. But from a low-level, it seems possible that a mismatch in the string reference ids might be a sign that they are effectively different..


texts:
Does anybody have a tutorial on how to download NASA NEX-GDDP downscaled GCMs dataset from http://opennex.planetos.com/ website as csv?
I´m not familiar with the use of Docker, containers, etc.
Any help would be thankful.


texts: Drop-down not showing on attribute table after I set a GDB domain?

texts: I created a file GDB with various feature classes that already have data in them.  I then set up a domain to the GDB, and when I open the attribute tables of each of my feature classes there is no-drop down to choose from the coded domain. I also tried adding a new field within the attribute table and there is still no coded domain.  Is this because my feature classes were pre-populated with data?


texts: Please assist me on how to use ArcCatalog to confirm the spatial reference is correct and adjust if necessary.


texts:
I'm using OSX 10.6 and just downloaded and installed PostgreSQL 9.0.2-1 and pgRouting 1.0.5-1 for Postgres 9.0 from here.

I'm a confused what my next steps are though. I'd like to try loading a simple network into pgRouting and running a few different network routing analyses on it to test it out. However, I don't know where to begin.

How do I open PostgreSQL or pgRouting? I'm fairly comfortable with terminal commands and scripting, but as I installed precompiled versions, I don't know where anything is or what to get at. I have looked at the pgRouting documentation, but I didn't see anything similar to my question.


texts:
I have been searching for similar problems however i did not find any that could fix my error.
i have GRASSGIS 7.0 running on MACOSX El Capitan. No problems at all until i tried to run i.pansharpen on Landsat 8 images.

this is how i set it (no optional values)

it tooks a very long time almost freezing on

  i.pansharpen --overwrite red=LC80340432016029LGN00_B3@PERMANENT
  green=LC80340432016029LGN00_B2@PERMANENT
  blue=LC80340432016029LGN00_B1@PERMANENT
  pan=LC80340432016029LGN00_B8@PERMANENT output=2016_pan_ method=brovey
  Performing pan sharpening with hi res pan image: 15.000000
  Histogram matching...
  Histogram matching...
  Histogram matching...
  Calculating Brovey transformation...

then after a while it goes like this:

  Performing pan sharpening with hi res pan image: 15.000000
  Histogram matching...
  Histogram matching...
  Histogram matching...
  Calculating Brovey transformation...
  Assigning grey equalized color tables to
  output images... Traceback (most recent call last):   File
  "/Applications/GRASS-7.0.app/Contents/MacOS/scripts/i .pansharpen",
  line 452, in
      main()   File "/Applications/GRASS-7.0.app/Contents/MacOS/scripts/i .pansharpen",
  line 319, in main
      flags="e", color='grey')   File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/pytho
  n/grass/script/core.py", line 395, in run_command
      return handle_errors(returncode, returncode, args, kwargs)   File "/Applications/GRASS-7.0.app/Contents/MacOS/etc/pytho
  n/grass/script/core.py", line 313, in handle_errors
      returncode=returncode) grass.exceptions.CalledModuleError: Module run None ['r.colors', '--q', '-e', 'color=grey', 'map=2016_pan__red']
  ended with error Process ended with non-zero return code -11. See
  errors in the (error) output.

after that the output is created empty...

i do not understand what is the error, and how i could fix it, i have been extensively searching for a similar problem but couldn't find anything.


texts: migration of irregular polygons to grids/blocks

texts:
How do you convert irregular polygons to grids (blocks) of regular/square shapes. I am trying to convert irregular shaped mining concessions to regular blocks in ArcGIS

Am using arcgis 10, arcinfo
I have created square grids of 45seconds for the whole country. Now I want to make sure all irregularly shaped mining concessions are converted to match the grids. Obviously a concession will lose/gain some areas.


texts:
I have a list of tweets on a torque map made with CartoDB. I want to add an additional layer that shows the sunlight around the world at certain times, similar to the Twitter #sunrise map found here.

Does anyone know how to add a layer like this?


texts: Performing in-memory difference operation between vector layers in PyQGIS

texts:
I've been working on a project using QGIS's Python API.  I need to create an in-memory layer that is defined as the difference between two polygonal vector layers, but I'm having difficulty figuring out how to do this for a standalone script.

I've looked at the processing plugin, but it seems to require a file output.  I've seen that some people use the runandload() function (How to load memory output from QGIS processing?), but this seems to be specifically for the QGIS GUI console, and I'm needing to generate the output layer in memory.


texts: Research papers on SWE (Sensor Web Enablement) and SOS (Sensor Observation Service) with security

texts:
I am looking for research papers (conference/journals) on the topic of SWE (Sensor Web Enablement) and SOS (Sensor Observation Service) also dealing/including security aspects.

I have already found some sources (after cleaning up, will post them here too) but you could post any expert material or websources where to search.


texts: Check whether values are unique in a column or not

texts:
How can I check if a value in a cell of a specific column is unique in the respective column or not.

e.g. I have a column called OBJECTID. Some of the ID's are unique whereas others exist twice.

I would like to create a new column which indicates whether a value is unique or not.


texts: Anyway to automatically cache the entire basemap in ArcMap?

texts:
At my job we are developing Add-Ins for ArcMap. We are expecting to do some demos for clients in the near future and to make the demos go as smooth as possible we would like to build up the local cache for the various basemaps found in the "Add Basemap" window.

Right now, it seems that the only way to build up the cache is for me to open ArcMap and start manually scrolling around the dataframe so that it downloads chunks of the basemap.

Is there anyway I could quickly get the entire basemap cached for every section at every scale level?


texts: NDVI image contains 500 cells or more using ArcMap

texts:
I want to make a map that shows my NDVI image (with 6 classess, natural breaks that I reclassified already)) that only contains zones with 500 cells or more and that are color coded to the area of each zones.

How would I do that? using which tools or function in arcmap? Im using ArcGIS 10.6

arcpy will help too


texts: ClientProperty of a top element not showing

texts:
On app-schema configuration file I have several elements with ClientProperty which shows the property on a WFS GetFeature response but when a ClientProperty is defined for a top element like the SIGMET defined below the id is returned but the two ClientProperties reportstatus and permissibleusage are not! Any ideas why ClientProperty doesn't work on a top element?

    <FeatureTypeMapping>
        <sourceDataStore>iwxxm_sigmet_postgres</sourceDataStore>
        <!-- sourceType is a database table or view -->
        <sourceType>sigmet</sourceType>
        <targetElement>iwxxm:SIGMET</targetElement>
        <attributeMappings>
            <AttributeMapping>
                <targetAttribute>iwxxm:SIGMET</targetAttribute>
                <idExpression>
                    <OCQL>gmlid</OCQL>
                </idExpression>
                <ClientProperty>
                    <name>reportStatus</name>
                    <value>reportstatus</value>
                </ClientProperty>
                <ClientProperty>
                    <name>permissibleUsage</name>
                    <value>if_then_else(isNull(permissibleusage), Expression.NIL, permissibleusage)</value>
                </ClientProperty>
            </AttributeMapping>

I'm getting back
<iwxxm:SIGMET gml:id="a946ee46-2cfb-4e65-8c07-e7f5f8670af01">

but it should be something like
<iwxxm:SIGMET gml:id="a946ee46-2cfb-4e65-8c07-e7f5f8670af01" reportStatus="NORMAL"  permissibleUsage="OPERATIONAL">


texts: GeoServer Image Mosaic - format date/time for TimestampFileNameExtractorSPI in indexer.xml

texts:
I'm trying to build mosaics from NetCDF files using the indexer.xml file instead of the *.properties files so we can take advantage of the NetCDF auxiliary file. (Using GeoServer 2.14.2 btw)

My typical filenames look like: filename_201903031503.nc

And using the old properties files my timeregex.properties is:

regex=[0-9]{12},format=yyyyMMddhhmm

However I can't figure out how to pass the format parameter when using the indexer.xml file, I started with the following collector configuration:

<collectors>
  <collector name="timeFrameCollector">
    <value>format=yyyyMMddHHmm,regex=[0-9]{12}</value>
    <spi>TimestampFileNameExtractorSPI</spi>
    <mapped>time</mapped>
  </collector>
</collectors>

As suggested on the following page: http://osgeo-org.1560.x6.nabble.com/ImageMosaic-indexer-xml-time-domain-custom-dimensions-td5144486.html

Unfortunately that gives me a Failed to create reader from file:data/supplier/product error.

As a temporary solution I can rename the files to this: filename_20190303T150300Z.nc

Then the following XML works:

<collectors>
  <collector name="timeFrameCollector">
    <value>[0-9]{8}T[0-9]{6}Z</value>
    <spi>TimestampFileNameExtractorSPI</spi>
    <mapped>time</mapped>
  </collector>
</collectors>

I'd really like to avoid renaming the files, they are supplied to us and I'd prefer to do zero/minimal processing. Is there a way to pass the date/time format to the collector using the indexer.xml file?


texts: Displaying WMTS-Layer in custom projection (EPSG:25832) using OpenLayers?

texts:
I have the URL of a WMTS server which provides map tiles in EPSG:25832 (ETRS89 / UTM zone 32N; I'm sorry I cannot post it as it contains user credentials).

If I connect to this server with QGIS I can select a tileset (ADV_25832_Quad) and it is displayed properly. The image below (left) shows the result with a vector layer (red border) for verification.

Now I tried to use this WMTS server in OpenLayers by combining these two examples:

WMTS example: https://openlayers.org/en/latest/examples/wmts.html
Custom projection example: https://maps.omniscale.com/de/examples/openlayers3_25832

In OpenLayers the WMTS layer is not displayed properly (see image below on the right).

I think is has to to with the projection extent or the resolution of the tilegrid (I got the projection extent from https://epsg.io/25832).

But how does QGIS get these values and how can I get them for OpenLayers?

Here's my code:

var projectionName = 'EPSG:25832';
proj4.defs(projectionName, '+proj=utm +zone=32 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs');
var projection = ol.proj.get(projectionName);

var projectionExtent = [ -1877994.66, 3932281.56, 836715.13, 9440581.95 ];
var size = ol.extent.getWidth(projectionExtent) / 256;
var resolutions = new Array(14);
var matrixIds = new Array(14);
for (var z = 0; z < 14; ++z) {
    resolutions[z] = size / Math.pow(2, z);
    matrixIds[z] = "ADV_25832_Quad:" + z;
}

var map = new ol.Map({
    layers: [
        new ol.layer.Tile({
            opacity: 1,
            source: new ol.source.WMTS({
                attributions: '...',
                url: 'https://...',
                layer: "DOP_20_C",
                "matrixSet": "ADV_25832_Quad",
                projection: projection,
                tileGrid: new ol.tilegrid.WMTS({
                    origin: ol.extent.getTopLeft(projectionExtent),
                    resolutions: resolutions,
                    matrixIds: matrixIds
                }),
                style: 'default',
                wrapX: true,
            })
        }),

        new ol.layer.Vector({
            source: new ol.source.Vector({
                url: 'bw25832.txt',
                format: new ol.format.GeoJSON()
            }),
            style: new ol.style.Style({
                stroke: new ol.style.Stroke({
                    color: 'red',
                    width: 2
                })
            })
        })

        ],
        target: 'map',
        view: new ol.View({
            projection: projection,
            center: [506102,5384251],
            zoom: 8
        })
      });

EDIT

Here are the top 2 levels of the tile matrix from the getCapabilities response:

    <TileMatrixSet>
              <ows:Identifier>ADV_25832_Quad</ows:Identifier>
              <ows:SupportedCRS>urn:ogc:def:crs:EPSG::25832</ows:SupportedCRS>
              <TileMatrix>
                    <ows:Identifier>ADV_25832_Quad:0</ows:Identifier>
                    <ScaleDenominator>1.747132075089743E7</ScaleDenominator>
                    <TopLeftCorner>-46133.17 6301219.54</TopLeftCorner>
                    <TileWidth>256</TileWidth>
                    <TileHeight>256</TileHeight>
                    <MatrixWidth>1</MatrixWidth>
                    <MatrixHeight>1</MatrixHeight>
              </TileMatrix>
              <TileMatrix>
                    <ows:Identifier>ADV_25832_Quad:1</ows:Identifier>
                    <ScaleDenominator>8735660.375448715</ScaleDenominator>
                    <TopLeftCorner>-46133.17 6301219.54</TopLeftCorner>
                    <TileWidth>256</TileWidth>
                    <TileHeight>256</TileHeight>
                    <MatrixWidth>2</MatrixWidth>
                    <MatrixHeight>2</MatrixHeight>
              </TileMatrix>


texts:
I've tried setting POSTGIS_GDAL_ENABLED_DRIVERS=ENABLE_ALL as an environment variable. As well as running

SET postgis.gdal_enabled_drivers = 'ENABLE_ALL'; from psql.
-- (With a pg_ctl restart after each)

I'm still getting 0 rows when I query

select short_name, long_name from st_gdaldrivers();

Postgres: PostgreSQL 9.4.1 on x86_64-apple-darwin14.1.0, compiled by Apple LLVM version 6.0 (clang-600.0.56) (based on LLVM 3.5svn), 64-bit

PostGIS: POSTGIS="2.1.7 r13414" GEOS="3.4.2-CAPI-1.8.2 r3921" PROJ="Rel. 4.8.0, 6 March 2012" GDAL="GDAL 1.11.1, released 2014/09/24" LIBXML="2.9.2" LIBJSON="UNKNOWN" RASTER

GDAL: (from gdalinfo --version) GDAL 1.11.1, released 2014/09/24

Is there something I didn't install?


texts: Remove duplicate values in an attribute using regex and fme

texts:
I have an FME process that is duplicating some values in an attribute. When I try to use regex in FME to clean it up it is deleting the separator (see screengrab below). Any ideas on how I can change the regex or process in FME to just get the required values.
So from
164511|119958|164512|164511|119958|164512|164511|119958|164512
161634|161635|161636|161634|161635|161636|161634|161635|161636|161634|161635|161636|161634|161635|161636|161634|161635|161636
1114|1114|1114|1114|1114

I just want
164511|119958|164512
161634|161635|161636
1114

The value in AttributeManager from FME is @ReplaceRegEx(@Value(Name_Ids),"(?:^|\G)(\b\w+\b),?(?=.*\1)\|",\1)
I have tried to replace \b),?(? with \b)\|(?to use the | seperator but that doesn't work either. This is based on https://superuser.com/a/1343524/639219
The regex tester is https://regex101.com/r/1YhC59/1


texts: Is the GEO tag for QR codes lat/long or long/lat

texts:
I just found out that you can add a geographic component to a QR code.

I also just found you can create one on the Google Chart API using this:

https://chart.googleapis.com/chart?cht=qr&chs=400x400&chl=geo:30.3,-97.6

That creates me a QR code that looks like this:

So far, so good. And when I read that back on my iPhone using the QR Code Reader app, it sends me to a Google map in the correct location (Austin, TX).

The thing is, the string used to create this uses lat/long. If I wrote the same data to a KML file it would expect long/lat.

So the question is, what spec do I use for QR codes:

Is the QR spec actually lat/long?
Is the Google Chart API getting it mixed up?
Is the QR reader app getting it wrong?

My money is on the QR reader messing up - but I'm not sure. Does anyone have a definitive answer?


texts:
I have several synoptic measurements (every few months) of a gravity field I would like to interpolate by kriging. Each group of measurements comprises about 45 measurement locations, and spatial correlation is relatively poor. I would like to use a single semivariogram model for all of the interpolated surfaces, rather than fitting a unique model for each measurement time. I think this is appropriate because the changes over time are caused by the same underlying phenomena (groundwater level change), and the spatial correlation should be unchanging with time.

Is there any problem with "stacking" the semivariogram data? I would have to calculate the semivariance-distance pairs independently for each measurement time, but could then combine all of them on a single semivarogram plot and fit my model to all of the data pairs. I'm primarly an arc user, but I think I could create the semivariogram model outside of arc (e.g., MATLAB) and simply specify the nugget/sill/major range for spatial analyst. Thanks-


texts: Unable to add Time Attribute to Network Dataset in ArcGIS Pro

texts:
I'm using ArcGIS Pro 2.

I want to setup service areas but I'm unable to use the time attribute. I want to use the cutoff in minutes to create service areas but the units for time say it's unavailable.

I have a network dataset. I have 'Time' as part of my attribute table. I don't know how to use this.


texts: GeoPandas cumulative value sum for intersecting polygons

texts:
I have a geojson file for all of the tax areas in the state of Utah. I read it into a GeoPandas dataframe. Now each row represents a tax area polygon.

Many of the tax area polygons overlap each other. For example, a county polygon serves as the foundation layer, and then a city polygon will be placed on top (the county polygon contains city polygons).

Each tax area polygon has a tax rate field. My objective is to add up tax rates where these polygons intersect, so in the end, one could hover over an area, and get the "top most" tax area's name, but have a cumulative sum of the intersecting tax area's rates.

I've looked at:

Projecting attribute from one polygon layer to another in QGIS?
Calculating area weighted average of attribute in one layer by shape in another using ArcGIS for Desktop? and
Finding all polygons overlapping polygons from another layer using MapInfo?

but these use other GIS software, and this is my first time using GeoPandas.

My ultimate goal is to output a new geojson file that has each tax area polygon but with the cumulative rate instead of each individual rate.

Any advice on how to accomplish this?


texts:
When I add inn an iterator into my model, the first tool (which in this case is "make XY event layer") keeps running in a loop. This occurs even if the iterator is not even connected to anything. Once I remove the iterator it seems to work as normal.

Any idea why?


texts: I have a university assignment that requires me to find an optimum location for a land fill site. There are various parameters I need to map. I have managed to map all the ones within a distance of something (using buffering and overlay techniques as in the image), but I have no clue how to define an area outside of a point. I need to define areas 20 miles outside of any major urban area, 20 miles from a major water feature and 10 miles from any national parks. I have all the base information (urban area, water features, national parks) in layers. I have a full licence with all additional add on's at my disposal - I just don;t know which to use or how]1


texts: Writing trigger which calculates driving distance?

texts:
Trigger fires but doesn't update my table - COUNTY_VR_GEO_ROUTING.

Here's the output from pgadmin after updating a record -

UPDATE 2
Query returned successfully in 2 secs.

Can anyone give me thoughts on what may be going on?

PostgreSQL 9.6; windows os.

CREATE OR REPLACE FUNCTION vr_routes.usp_etl_gis_pgr_dist_tigger_fn()
RETURNS TRIGGER

AS $BODY$
BEGIN
      UPDATE "VOTER_ROUTING".vr_routes."COUNTY_VR_GEO_ROUTING"
      SET dist = agg_cost * 0.000621371 FROM (SELECT agg_cost FROM pgr_dijkstraCost(
           'SELECT gid as id, source, target, length_m AS cost FROM osmways',
            -- source
        (SELECT id FROM osmways_vertices_pgr
            ORDER BY ST_Transform(the_geom,2264) <-> ST_SetSRID(ST_Point(NEW.x1,NEW.y1),2264) LIMIT 1),
        -- target
        (SELECT id FROM osmways_vertices_pgr
            ORDER BY ST_Transform(the_geom,2264) <-> ST_SetSRID(ST_Point(NEW.x2,NEW.y2),2264) LIMIT 1),
        false)
         ) AS dijkstra; -- undirected

 RETURN NEW;

END;
$BODY$
LANGUAGE plpgsql

Here's my trigger

CREATE TRIGGER usp_etl_gis_trigger_county_vr_geo
AFTER INSERT OR UPDATE
ON "VOTER_ROUTING".vr_routes."COUNTY_VR_GEO_ROUTING"
FOR EACH ROW
EXECUTE PROCEDURE vr_routes.usp_etl_gis_pgr_dist_tigger_fn();


texts: Is there a source for leaf-off aerial imagery?

texts: I am looking for leaf-off imagery for Oklahoma, Kansas and Texas.  1m leaf-off imagery would be ideal, although any substitute would be helpful.  I need fairly high resolution imagery to distinguish conifer from deciduous trees--Landsat would likely be too coarse.  How can I obtain leaf-off imagery for KS, OK, TX?


texts: Using GeoServer CSS styling, can you switch off transformation at a certain zoom scale?

texts:
I have came across similar questions and read through the documentations, so sorry if I missed something. I'm looking for a way to turn of transformations based on a zoom scale using CSS. For an example, I have written the following code to cluster points together using the point stacker:
[@sd >= 100k]{

  transform: vec:PointStacker(cellSize:30, outputBBOX:env('wms_bbox'), outputWidth:env('wms_width'), outputHeight:env('wms_height'));
  mark: symbol("circle");
  :mark{

        fill:  [categorize(
         count,
         1,
         1, '#4866A3',
         25, '#26ADF2',
         50, '#53F08E',
         100, '#21E86D',
         200, '#FF9D35',
         400, '#EA4638')];

        stroke: black;
        size:  [categorize(
         count,
         1,
          1, 3,
         10, 5,
         25, 10,
         50, 15,
         100, 20,
         200, 30,
         400, 50)];

}
}

After a user has zoomed in (@sd < 100k), I would like to show single points categorized/recoded on a field, for an example:
[@sd < 100k]{

   mark: symbol(circle);
   :mark {
      fill: recode([HomeAddress_Confidence], 1,'#1A62C9', 2,'#F0FBF7', 3,'#0398EE', 4,'#3EEC83', 5,'#FFBE73', 6,'#D35230', 7,'#C90000') ;
      size: 5;
   }

}

It seems that although PointStacker are within one rule it takes over the whole css and can't use attribute in other rules. I also included the different modes using @mode 'Eclusive', 'Flat', 'Simple' and 'Auto'. A possible work around is to create 2 separate CSS styles and create a group layer, but would like to avoid that if possible.


texts: Seeking Seabed Survey Data Model (SSDM)?

texts:
Many of the offshore oil/gas companies are adopting an initiative from Shell for delivery of spatial data.  This is called an SSDM, which is really an ESRI file geodatabase with a standardised structure of datasets, and featureclasses.

An empty template (geodatabase) in the SSDM structure is freely available from the OGP website, but it would be really good if anyone knows of an example geodatabase which is fully verified so we can use it as a 'how to' reference guide.

Does anyone know if such a thing exists?


texts:
I know there is a way to download a local copy of OSM. I would like to be able to create a second instance of that copy, and to modify the second instance (adding planned roads, etc). Quick Googling suggests that I can either use ArcGIS to make edits.
Is there a way to make MASS edits (in a geoprocessing sense) to an OSM dataset, ie by spatial joining count data to links?
Will I need to do something like converting the OSM format file to a shapefile, make edits, and convert back?


texts: QGIS GUI: How to get independant windows?

texts:
I am working with the georeferencer plugin, and I would like to have it working in fully independent windows.

Is it possible ?

Thanks for your help.

Context:

I am working with qgis 2.4 under Ubuntu 12.04.

I spot a dock/undock button on the top right corner of the plugin window.
However even when undock is pressed, the windows remains linked to the qgis window and cannot be moved freely or displayed in full screen.

Here is as overview of the windows layout


texts: How to retrieve bus line including stop names?

texts:
I don't know much about OSM, OverpassTurbo, but would like to be able to query OSM about a given bus line, and retrieve the line and all the station names.

Then, I'll export the data as GPX so I can use them on my smartphone with no Internet access.

The following works, but I don't know how to customize it so it also retrieves the station names:

(
  relation["network"="NameOfBusLine"]["ref"="123"]({{bbox}});
);


texts: Some companies receiving GPS data from vehicles when GPS cannot send any data?

texts:
I know that some companies offer tracking services for vehicles such as live GPS tracking they can access any time, but as you can't use the GPS protocol to send any data i wonder how the do it... I thought it can be archived with some attached GPRS functions to the GPS but the cellphone network is less reliable than GPS for some remote areas.

So how that's possible?

Is there any additional network attached?


texts:
I want to give feedback to the developer of a non-core QGIS plugin.

(Numerical Vertex Edit by Cedric Moeri, who seems to have no webpage for his plugins)

I joined OSGeo Wiki but hardly know my way around it, or whether that is a good place to start.


texts: Get the attribute of the nearest line to a polygon

texts:
In QGIS, how could I get the attribute of the nearest line to a polygon with an expression (or virtual layer)?
I wish to create a new field in the polygon layer(in this case the parcels) with the names of the streets (the line layer) that are nearest of them. I thought I could use aggregate but I didn't find any nearest function to go with it.


texts: Spatial Queries in Leaflet, select polygons from (A1,A2...An) which intersect with polygon (B)

texts:
I am looking for an online solution for leaflet where I can select and export (console log, geojson, addLayer) a group of polygon from a geojson file which intersect with a drawn polygon (with leaflet draw). Please see the picture:

I tested Esri leaflet but for some reason it is not working form me. I am not sure if it requires an online esri license or I did something wrong (my layer was a geojson not an esri source layer), however I did not get any error messages and also nothing happened. Could anybody confirm this plugin works free without the necessity of any esri related license or data?. If not, my second question is can this question be solved with turf javascript? I know there is turf intersect: var intersection = turf.intersect(poly1, poly2); but this is between 2 polygons and also the end result is new geometry and not a selection of existing ones.

I also find this link which looks useful, but doesn't work for me and I guess the end result would be a new geometries of the intersection:
Intersecting FeatureCollections in turf.js?


texts:
I'm very new to the GIS world, so my problem could be very stupid, but I'm going to try anyway.

Objective

Given a list of coordinates in longitude/latitude, stored as a .csv file, I want to create a buffer with distance of x km around those coordinates. (In case this is relevant, these coordinates are locations of the societies in the SCCS dataset. These societies scatter all around the world.)

Problem

I'm only able to create buffers with distance in radial degrees, but I want to do them in (kilo)meters. I'm aware of the following Q&A's, but following what's suggested therein doesn't seem to solve my problem:

Understanding QGIS buffer tool units?
What's the Workflow for buffering in meters/km using WGS84 layers?
How to create a polygon (fixed width buffer) around a polygon in QGIS?
How to change the unit of measure from degrees to meters in QGIS?

What I have tried...

I start QGIS (v2.4), and Add delimited text layer
choose my .csv file and the x/y fields, click OK
select WGS 84 as my CRS (also tried NAD83 but didn't make a difference), click OK
right-click on layer and select Save As..., in the dialog box, choose the following and click OK:

add sccs_meter.shp as layer to the current project
Select Vector > Geoprocessing Tools > Buffers > set buffer distance as 10

The result is a buffer with radius of approximately 1068km, which suggests that buffer distance is 10 degrees rather than 10 meters.

What did I do wrong here?


texts: ST_Equals not matching with exact geometry

texts:
I followed many question posted on the same subject but unable to solve my problem.

I have same kind of tables in 2 different database (tables with linestring geometry), I need to update id from one table to another when it finds equal geometry in another table.

table1 from DB1 looks like bellow.

     Column      |           Type           |
-----------------+--------------------------+
 old_id          | bigint                   |
 road_id         | bigint                   |
 road_name       | character varying(255)   |
 geo             | geometry(Geometry,4326)  |
Indexes:
    "roads_final_pkey" PRIMARY KEY, btree ("old_id")
    "roads_final_wkb_geometry_geom_idx" gist (geo)

and table2 of DB2 is as bellow.

   Column    |          Type           |
-------------+-------------------------+
 id          | bigint                  |
 object_id   | bigint                  |
 geoproperty | geometry(Geometry,4326) |
Indexes:
    "geofeature_pkey" PRIMARY KEY, btree (id)
    "geoproperty_idx" gist (geoproperty)
    "object_id_idx" btree (object_id)

Before updating I just checked how many of them are matching, like bellow.

select count(*) from (select * from dblink('port=5432 dbname=db2 user=usr password=pass','SELECT geoproperty from table2' as t(geoproperty Geometry)) r,table1 t1 where st_equals(t1.geo,r.geoproperty);

But it is returning only 10-15% of matches where I was expecting more than 90% match because they are same geometry copied to DB2 and updated object_id.

I also tried with st_snaptogrid,st_buffer and st_astext with geo columns, But they are taking huge amount of time compared to regular st_equals method.

But when I cross checked with non matched id's geometry from table1 with geometry in table2 with normal

select * from st_equals(st_geomfromtext('LINESTRING(.line from table1.)'),st_geomfromtext('LINESTRING(.line from table2.)'))

It is returning true but not matching with first (count/select) query.

Is there anything I am missing ?


texts: Performing reverse-clip/difference() point GeoJSON from polygon GeoJSON in QGIS?

texts:
I am trying to clips the points outside of a polygon GeoJSON.

In QGIS(2.18.2) vector-> geoprocessing tools -> difference deletes the geometry but keeps the attributes of the points.

The hacky way is to use vector -> research tools -> select by location to select all points inside the polygons and then delete them.