texts:
I'm having trouble exporting my map in TileMill (v0.10.1). I'm using OSM Bright with a PostGIS database connection to import OpenStreetMap data. When exporting an SVG or PDF through the user interface in TileMill, the export fails more often than not with a SIGABRT error. The console shows the following error:

[tilemill] [tilemill-ui] Rendering file
[tilemill] [tilemill-ui] libc++abi.dylib: terminating
[tilemill] Export process failed with signal 'SIGABRT'

I then tried exporting through the command line by running ./index.js export OSMBright ~/Documents/MapBox/export/cityhall2.svg --format=svg --width=2000 --height=2300 --bbox="-79.8786965, 43.2604005, -79.8702965, 43.2516005", with a nearly identical error.

It usually exports a blank file, although once in a while it does export nearly what I want (with zoom levels and symbology slightly off, but still usable).

Has anyone had this problem before? Thanks for the help!


texts: How to convert float raster to integer raster?

texts:
I am using Normalize difference vegetation index (NDVI) derived from the landsat satellite band 3 and band 4 that have float values ranges from -1.0 to 1.0. I want to see the attribute table that is not possible in QGIS.
I tried following expression in raster calculator of qgis but it did not work:

Int("float_raster" * 100)

the reason to multiply with 100 is to preserved the 2 digits at right level in attribute table that can be seen by again dividing by 100 in attribute table.

How can i scale see the better contrast in colour palette using qgis 2.40


texts: Complicated vector layer's style in OpenLayer 2

texts: Is it possible to style vector layer so there would be picture and label with its own background color for a feature? Something like this:


texts: Identifying "Must not connect point more than one line feature" errors in electrical network?

texts:
I have attached an example image to show the problem: This is electrical network data single service point to connect more than one Service pole. Is it available any topology rules?
I would like some suggestions regarding this error .


texts: How to display Tile Coordinate or Quadkey in QGIS map?

texts:
I am trying to create the following type of map with QGIS:

For the zoom level displayed, I want to show the contour of the tiles and the tile coordinates.
What is a quick and simple way to do it ?
QGIS: 3.10


texts: How can I convert point features into a raster that the original point will be the center of the new cell in ArcGIS?

texts: I have a problem in ArcGIS, So, I converted rasters to points. After that, I ran extract values to ponts function. Now I want to convert my point features back to raster with the original resolution. The problem is new cells are not at the same location as the old ones because I ran point to raster function. When ArcGIS convert ponts to rasters the points are not in the central of the new cells. Just how you can see on picture.


texts: Converting polylines to raster, assigning value of total length of lines within cell?

texts:
I have a polyline shapefile representing a road network. I want to rasterise this data, with the resulting values in the raster showing the total length of lines that fall within the raster cell, similar to this question: Convert line shapefile to raster, value=total length of lines within cell - however I wish to use ArcGIS 10.2. The data is projected so the units will be metres.

I have attempted to do this via the polyline to raster tool (however the value assigned is the length of the whole segment) and also some spatial analyst tools such as line density (which requires a radius). I'm sure I am missing an easy way to do it but I can't seem to find it.


texts: Mobile app for workflow/timesheets + GIS Asset collection required

texts:
We're looking for a solution to move our field operatives from paper based input to electronic field input on tablets to record workflow and map assets. Work involves land drainage (working on watercourses - reforming/reprofiling, desilting, flailmowing, improvement works etc.)

We need to make it as simple to use a possible so we want an integrated form + mapping solution. We've looked at things like ESRI collector and Survey123 but they're too basic for our needs.

Example workflow for an operative:

Start of day, gets job on tablet.

Use tablet to log start of work, vehicle start mileage etc.

Complete mandatory vehicle checks, risk assessments etc. (need to be recorded on tablet).

Does job, vehicle tracked on GIS map, throughout job operative can input asset data on to GIS.

Ends job, time recorded on tablet + additional info.

This is a very simplified example and more info will need to be captured throughout. Ultimately all the info needs to get back to head office to be verified and processed.

Big picture is that we want to move away from paper timesheets etc. as the app should be able to provide the info to generate this and asset info collected in the field will go on to the central GIS.

Unfortunately don't really have any in house capability to develop this ourselves. We do have good desktop GIS experience and knowledge. We're UK based and hoping to be able to outsource the app development. Any suggestions welcome?


texts:
Could anyone just clear up a question regarding openstreetmap data. Are there different OSM data sets for the various set scales, or do the data get generalized on-the-fly when the map is zoomed out to a smaller scale? I understand that the many objects are hidden in smaller scales, but are there also different generalized datasets for the smaller scales as well?

Thanks,

Robert


texts: Create Custom Tool with ArcObjects SDK and Visual Studio 2013/2015

texts:
I'm researching building a couple of custom tools.  I downloaded and installed ArcObjects SDK 10.4 to integrate it into Visual Studio 2013 and 2015, but I'm having trouble finding any good examples on creating the tools.  Our client wants one tool that will search for a value from a specific field("Loc_Num") within a specific set of layers and then highlight the feature it finds.  They also want another one to do an address lookup, but apparently the built-in tool is not what they are wanting.

Can anyone could point me to a good reference on building these within VS 2013 or 2015?


texts: What sources of granular wind data over the continental US are available?

texts:
I am looking for publicly or commercially available wind data at altitudes ranging from 0 feet to 40,000 feet over the continental United States.

Which data sources should I look at? Which one is available at the finest resolution (both spatially and temporally) possible? How was this data measured or modeled?


texts: QGIS print composer and label with svg background error

texts:
I have a problem with labels in QGIS. I'm using a background svg for some labels (street numbers). In QGIS, it's all ok but when I use the composer there is a bug (I suppose), because I see some strange artifacts.

image

svg

If I disable the layer with the street numbers, or if I remove the svg background, the artifacts disappear. The problem is reproducible only with svg, if I use the rectancle or another option, I don't see any problems.
The artifacts can be seen in the composer preview too.

Can someone confirm this problem?

Tested with QGIS 2.6.1 and 2.8.1 Nightly


texts: QGIS "effects" not displaying in print composer

texts:
I'm having an issue where effects created under Layer Properties -> Style -> Draw Effects/Customise Effects do not display in the print composer window.

See screenshot below - in the main QGIS window I have an outer glow effect applied to the workings svg symbol. However, once a print composer is created, with a map displaying that area, the effect disappears, just leaving the standard black svg icon.

I'm sure I've previously had maps in the print composer honouring custom effects, but I may be mistaken.

Any advice?

I'm using QGIS 2.18.5.


texts:
I am using the following sequence of commands to extract some data from open street map to import to Postgres:

osmconvert -b="-7.57216793459,49.959999905,1.68153079591,58.6350001085" -o=/data/UK.poi.osm /data/europe-latest.osm.pbf

osmfilter /data/UK.poi.osm \
    --keep="highway=bus_stop amenity=clinic =dentist =hospital =doctors =place_of_worship =bar =pub =restaurant =cafe =school =university =pharmacy =cinema =nightclub shop=convenience =supermarket =parking_space =parking =retailaddr building=hotel railway=station addr:housenumber= addr:housename= addr:flats= addr:street= addr:place= addr:city= addr:postcode= url= phone= description=" \
    --drop-relations --drop-author --drop-version > /data/UK_filtered.poi.osm

osmconvert --all-to-nodes --csv="@id @lon @lat amenity name building highway railway url phone description addr:housenumber addr:housename addr:flats addr:street addr:place addr:city addr:postcode" --csv-headline --csv-separator=',' -o=/data/UK.poi.csv /data/UK_filtered.poi.osm

What I get back is a list of empty rows like:

5710782475,-6.0900342,54.4808682,,,,,,,,,,,,,,,

and also misformed lines like:

1000000100018169,-4.6536670,54.0737920,,Radcliffes,yes,,,,+44 1624 822271,"Craft butcher  deli & fishmonger,6,,,Malew Street,,,

where there is a missing double quotes at the end, invalidating the csv.

How can I fix these issues?


texts: Rule based coloured background for labeling in QGIS

texts:
I am a newbie working on QGIS 3.6 / Windows 10. I need to label a polygon layer which is categorised by colour.
(pls refer to attached pic as a graphical reference)

The labels need to have a background which is also categorised by colour (red, orange, yellow and white) according to values in the attribute table of the shapefile.

My understanding is that the map as pictured has been created with an older version of ArcView...

Is it even possible to achieve something similar in QGIS?


texts: Appending polygon features from file geodatabase table to SDE table?

texts:
I am facing an issue while appending the feature class to ArcSDE feature class.

All the I am trying to append the gdb feature classs data, I can see that the geometry is either changed and some cases they are NULL. Here is the code which I tried.

    arcpy.CopyFeatures_management(inputTable, str.format("{0}\\{1}",GdbPath, "FeatureTempTable"))
    arcpy.MakeFeatureLayer_management (str.format("{0}\\{1}",GdbPath, "FeatureTempTable"), str.format("{0}\\{1}",GdbPath, "FeatureTempTableLayer"))

     arcpy.SelectLayerByAttribute_management(str.format("{0}\\{1}",GdbPath, "FeatureTempTableLayer"),  "NEW_SELECTION", "OPTYPE='I'")
    if int ((arcpy.GetCount_management(str.format("{0}\\{1}",GdbPath, "FeatureTempTableLayer"))).getOutput(0))>0:

        arcpy.CopyFeatures_management (str.format("{0}\\{1}",GdbPath, "FeatureTempTableLayer"), str.format("{0}\\{1}",tempGdbPath, "TempTable"))
      #In the Append below the shape is getting corrupted and sometimes becoming Null
        arcpy.Append_management(str.format("{0}\\{1}",GdbPath, "TempTable"), outputTable, "NO_TEST")

In Geodatabase I checked the shapes are good and even used the SHAPE@WKT in gdb. It seems good.
What I am trying is to select all the data in the feature class in the gdb where OPTYPE is I, are selected and appended to the SDE feature class in Oracle database.

What should I do?


texts:
I have make an interactive map with Leaflet and when I click in any point at the basemap  it displays an empty pop-up.

Here is The code
map.addEventListener('click', onMapClick);

    function onMapClick(e) {

    var latlngStr = '(' + e.latlng.lat.toFixed(3) + ', ' + e.latlng.lng.toFixed(3) + ')';

    //var BBOX = map.getBounds()._southWest.lng+","+map.getBounds()._southWest.lat+","+map.getBounds()._northEast.lng+","
    //+map.getBounds()._northEast.lat;
    var BBOX = map.getBounds().toBBoxString();
    var WIDTH= map.getSize().x;
    var HEIGHT = map.getSize().y;
    var X = map.layerPointToContainerPoint(e.layerPoint).x;
    var Y = map.layerPointToContainerPoint(e.layerPoint).y;
    var layerObj = {};
    layerObj["cite:Proj_Agadir"] = Projet;
    layerObj["cite:Projets_AssaZag"] =Projets_AssaZag;
    layerObj["cite:Proj_Inzegane"] = Projets_Inzegane;
    layerObj["cite:Proj_Taroudant"] =Projets_Taroudant;
    layerObj["cite:Proj_Tata"] = Projets_Tata;
    layerObj["cite:Essaouira"] =Projets_Essaouira;
    layerObj["cite:Export_Output3"] = Projets_Guelmim;
    layerObj["cite:Chtouka"] =Projet_Chtouka_Ait_Baha;
    layerObj["cite:SidiIfni"] = Projet_Sidi_Ifni;
    layerObj["cite:Tiznit"] =Projet_Tiznit;
    layerObj["cite:AEP"] = AEP;
    layerObj["cite:Education"] = Education;
    layerObj["cite:Piste"] =Piste;
    layerObj["cite:Sante"] = Santé;
    layerObj["cite:AGR"] =AGR;

    queryLayers = [];
    for (key in layerObj) {
    if(map.hasLayer(layerObj[key])) {
     queryLayers.push(key);
     }
    }
    var urlLayers = queryLayers.join();
    var URL = 'http://localhost:8086/geoserver/cite/wms?SERVICE=WMS&VERSION=1.1.0&REQUEST=GetFeatureInfo&LAYERS=cite:AEP,cite:Education,cite:Piste,cite:Sante,cite:AGR,cite:Proj_Agadir,cite:Projets_AssaZag,cite:Proj_Inzegane,cite:Proj_Taroudant,cite:Proj_Tata,cite:Essaouira,cite:Export_Output3,cite:Chtouka,cite:SidiIfni,cite:Tiznit&QUERY_LAYERS='+urlLayers+'&STYLES=&BBOX='+BBOX+'&FEATURE_COUNT=5&HEIGHT='+HEIGHT+'&WIDTH='+WIDTH+'&FORMAT=image%2Fpng&INFO_FORMAT=text%2fhtml&SRS=EPSG%3A4326&X='+X+'&Y='+Y;
    popup = new L.Popup({maxWidth: 1000});

    popup.setLatLng(e.latlng);

    popup.setContent("<iframe src='"+URL+"' width='1000' height='200' frameborder='0'></iframe>");

    if (urlLayers)

    map.openPopup(popup);

}


texts: Exporting thumbnails to jpg from ArcGIS feature class raster field with ArcPy?

texts:
I am trying to export rasters stored in a raster field in a gdb feature class to thumbnail-sized jpg with arcpy. The reason is so that I can attach them to the feature class to be displayed in a popup in ArcGIS online.

If I click on the raster field in ArcMap, the thumbnails are shown, but I am not sure how to export them with arcpy.

I have considered copying the raster to jpg, but what I really want is a down-sampled jpg (ideally less than 1mb each). I could use ExportToJPG, but then I would need to set up an mxd.

The @FelixIP suggestion to resample works well to create thumbnails. Here is a functional snipped to create thumbnails, and add the URL and thumbnail name to a table that could be used to show the thumbnail in AGOL.

#Export Raster to Thumbnail Tif

import os
import arcpy
from arcpy import env

rcat = "location of an unmanaged raster catalog"
pathsName = "name of output location of a table of raster paths"

#Export raster catalog paths to dbf table
arcpy.ExportRasterCatalogPaths_management(rCat, "ALL", pathsName)

# Add fields to store some useful fields to help manage thumbnail paths
arcpy.AddField_management(pathsName, "thumbURL", "TEXT", field_length=500)
arcpy.AddField_management(pathsName, "rasName", "TEXT", field_length=100)

#Use data access cursor to loop through table and create thumbnails
with arcpy.da.UpdateCursor(pathsName, seColumns) as cursor:
    for row in cursor:
        #Assign the raster path to a variable from the table
        rasPath = str(row[1])

        #assign input raster name to a variable
        inName = arcpy.Describe(rasPath).name

        #Create a URL for the thumbnail to be stored later
        tURL = "web server location" + inName + ".png"

        #Write thumbnail URL and raster name
        row[2] = tURL
        row[3] = inName
        #update the rows
        cursor.updateRow(row)

        #Create thumbnail to send to web
        outThumbLoc = 'tumbnail storage folder'

        #Create thumbnail name
        outThumb = os.path.join(outThumbLoc, inName + ".png")

        #Use resample to create thumbnails
        try:
            arcpy.Resample_management(rasPath, outThumb, tScale, "BILINEAR")
        except Exception as err:
            print rasPath
            print(err.args[0])


texts:
Here is the relevant portion of my mapfile:

CONNECTIONTYPE POSTGIS
  CONNECTION "dbname=contours host=localhost user=contour_xxxx password=xxx"
  PROCESSING "CLOSE_CONNECTION=DEFER"
  DATA "wkb_geometry from contours USING UNIQUE ogc_fid USING
contourinterval=40"

I want to limit the contours returned to only include contourinterval=40.
As you can see from the Postgres log (below), contourinterval isn't even
referenced in the query...

select
    "contourelevation",
    encode(ST_AsBinary(ST_Force2D("wkb_geometry"),'NDR'),'hex') as geom,
    "ogc_fid"
from contours
where wkb_geometry &&
ST_GeomFromText('POLYGON((-111.939697265625
33.4213636331582,-111.939697265625 33.4323499612832,-111.9287109375
33.4323499612832,-111.9287109375 33.4213636331582,-111.939697265625
33.4213636331582))',find_srid('','contours','wkb_geometry'))

I've tried all sorts of subselects and examples found on the web, nothing
works. BTW, I assume 'contourelevation' is pulled because it is my
LABELITEM. The data is USGS contour data.


texts: Performance issues with PostGIS ST_Intersects

texts:
I have a couple of tables in PostGIS that are big (~4 million rows). The geometry columns in both the tables are indexed. But I find a huge difference in query times between ST_Within and ST_Intersects.
Here is a stripped version of the tables:
Table Trips:

tripid
startgeom (point)
endgeom (point)

Table Paths:

tripid
geom(linestring)

These are the queries: (I have replaced the actual geometry value with a placeholder for privacy issues here)
Query 1:
explain analyse analyse select trips.tripid
from trips
where
st_within(trips.start_geom, '<geometry value>');

Bitmap Heap Scan on trips  (cost=154.00..18002.83 rows=1476 width=33) (actual time=38.553..1998.805 rows=3872 loops=1)
 Recheck Cond: ('<geometry value>'::geometry ~ start_geom)
 Filter: _st_contains('<geometry value>'::geometry, start_geom)
 Rows Removed by Filter: 1028
 Heap Blocks: exact=4846
 ->  Bitmap Index Scan on trips_start_geom_idx  (cost=0.00..153.63 rows=4429 width=0) (actual time=37.375..37.376 rows=4900 loops=1)
       Index Cond: ('<geometry value>'::geometry ~ start_geom)
Planning time: 0.122 ms
Execution time: 2002.607 ms

Query 2:
explain analyse select paths.tripid
from paths
where st_intersects('<geometry value>', paths.geom);
st_intersects('<geometry>', paths.geom);

Bitmap Heap Scan on paths  (cost=37.99..3099.98 rows=245 width=29) (actual time=10443.387..301686.132 rows=240341 loops=1)
 Recheck Cond: ('<geometry>'::geometry && geom)
 Filter: _st_intersects('<geometry>'::geometry, geom)
 Rows Removed by Filter: 1060363
 Heap Blocks: exact=628044
 ->  Bitmap Index Scan on paths_geom_idx  (cost=0.00..37.93 rows=735 width=0) (actual time=10178.236..10178.238 rows=1300704 loops=1)
       Index Cond: ('<geometry value>'::geometry && geom)
Planning time: 0.195 ms
Execution time: 301951.620 ms

The geometry used in both the queries is a small polygon. But there is a huge difference in execution times: ~2s to ~300s. Any idea why this is?
Update: Both tables have the exact same number of rows (~4 million).
Update 2: Query times after st_subdivide.
Query 3:
select paths_subdivided.tripid
from paths_subdivided
where
st_intersects('<geometry value>', paths_subdivided.geom);

Analysis before indexing:
Gather  (cost=1000.00..1698949.93 rows=504 width=28) (actual time=1.470..63191.646 rows=263736 loops=1)
  Workers Planned: 2
  Workers Launched: 2
  ->  Parallel Seq Scan on paths_subdivided  (cost=0.00..1697899.53 rows=210 width=28) (actual time=0.845..62986.318 rows=87912 loops=3)
        Filter: (('<geometry value>'::geometry && geom) AND _st_intersects('<geometry value>'::geometry, geom))
        Rows Removed by Filter: 1800743
Planning time: 1.519 ms
Execution time: 63393.985 ms

Analysis after indexing:
Bitmap Heap Scan on paths_subdivided  (cost=70.46..5578.92 rows=442 width=28) (actual time=624.927..311657.669 rows=263736 loops=1)
   Recheck Cond: ('<geometry value>'::geometry && geom)
   Rows Removed by Index Recheck: 1351624
   Filter: _st_intersects('<geometry value>'::geometry, geom)
   Rows Removed by Filter: 1005361
   Heap Blocks: exact=367315 lossy=332051
   ->  Bitmap Index Scan on paths_subdivided_geom_idx  (cost=0.00..70.35 rows=1325 width=0) (actual time=482.389..482.391 rows=1269097 loops=1)
         Index Cond: ('<geometry value>'::geometry && geom)
 Planning time: 0.404 ms
 Execution time: 311947.221 ms


texts: How to find good meter-based projection in PostGIS

texts:
I have a set of point geometries that I would like to cluster together using ST_ClusterDBSCAN.  Unfortunately, it doesn't seem to work for WGS84 and requires a meter-based projection.  Every tutorial that I can seem to find manually determines which projection is a good meter-based projection for their particular dataset.  Is there a generalizable way to do this automatically?

Basically what I'd like to do is take the convex hull of my points, take the centroid of that polygon and then query to find what's a good meter-based projection for that particular region of the world.  I'll then convert my points to this projection and run the clustering.


texts: I have a raster layer of which I would need to find the area (m2) per cell. Any suggestions for a simple approach to achieve this?


texts: Is it possible to make a multiple line coordinates in QGIS 2.12.1 map comsposer?

texts:
I'm trying to display a 2 line coordinates, where the thousands are on the top row, and the units are on the lower row.
I've tried to custom coordinates using the following expression:

left(@grid_number,3)  ||  '\n'  || right( @grid_number ,3)

But instead of a new line, I receive a font separating the 2 expressions that should be presented:


texts: How to customize and add additional features in GeoNode

texts:
I need some advice on how to start customizing my GeoNode installation. I tried reading up on GeoNode documentations on how to do it but it always boils down to customizing or adding to the Django framework. I am confused on which end should I edit first, should I focus on GeoNode only or should I read up on the Django framework?

Functions I'm trying to add:

attribute table querying
multiple feature highlighting
changing the built in map search function

I am still a bit new to this technology.


texts: Batch set value/calculate field inside Cursor

texts:
I am matching data using a cursor, where one of the final steps is filling in a new attribute, match, with a value calculated during the cursor. This matching is done to an ID field in the original data that is often shared by many points (i.e. there are 200 points with VID 8, and they will all end up with the same match value, 284). I am looking for a way to speed up my script by setting the match value on all common VIDs so instead of running the cursor on every point data, it is only run on the unique VIDs (unfortunately other steps in my script make dissolving this data into multi-point impossible). If this were not done inside a Cursor, I would select by VID attribute and then use the calculate field tool. Is there a way to batch set values or get around the schema lock issue for the calculate field tool inside a cursor?
Here is a sample of the code (simplified) if I could use Calculate Field:

rows = arcpy.SearchCursor("logbook") #Temporarily a SearchCursor, but I've tried with UpdateCursor as well.
for row in rows:
    if row.match > 0: #If the VID has already been matched, this will skip the calculations for it
        continue
    else:
        #(My matching and value calculating processes)#
        arcpy.SelectLayerByAttribute_management("logbook", "NEW_SELECTION", ' "VID" = ' + str(row.VID))
        arcpy.CalculateField_management("logbook", "match", int(match_val))
del row, rows

Any thoughts? Or am I doomed to have a script that takes much longer than necessary?


texts: How to normalize multiple PNG rasters with differing colormaps

texts:
I have a series of adjacent rasters which comprise a digital elevation model. Unfortunately they have been given to me as PNG files, rather than Esri Grids. The only key that I have to the actual elevation values is a separate legend PNG.

Each PNG file has an inbuilt colormap file, and they display correctly when drawn in ArcMap using the automated Colormap symbology option. The red line in this screenshot shows the boundary between two adjacent images:

However, Identifying within the dark green sections of each image shows that while the colors match (a,r,g,b = 255,0,225,0 on both the upper and lower images) the Color Index values do not match.

This means that if I try to Mosaic these images, there are discrepancies between the values. I'm using the Mosaic To New Raster tool with the Match Colormap option, but this seems to match the Color Index values rather than the colors(?).

Exporting the colormaps and comparing them in a text editor confirms that the Color Index vs a,r,g,b values do not match between the images, and that there isn't a simple shift between the values in each image:

Does anyone have any suggestions on how I can normalize (if that's the correct term) these adjacent images so they are using the same values? I can then mosaic and reclassify this single image to match with the legend.

(A cumbersome workaround would be to reclassify each individual image based on the legend, but I'd like to avoid this as it'll be hard to scale).

Thanks.


texts: QGIS how to patch a plugin on local machine

texts:
I want to make a minor modification to a plugin (to change a field width in Table Manager).

I edited the .py file and the .ui file and changed what I think are the right numbers but this seems to make no difference when I restart QGIS.

I observed that the .pyc file does not get up updated so I removed it but that does not make any difference.  Is the plugin cached somewhere?  I have been reading the docs for writing plugins and they seem to imply that changing the source and restarting QGIS is all that is necessary.

Of course I may have not understood the code and everything is working exactly as it should!


texts: Removing overlap in raster using Erase command line?

texts:
I am trying to erase a raster and is there a erase command for the Raster side?

I am using ArcGIS 9.3.1.

I have 2 rasters that are in the Landcover and there is overlap in between and would like to remove the overlap and then merge them later.


texts:
I recently received a collection of satellite images (GeoTiffs) and their indice. Originally, the indice came in a geographic projection (GCS WGS 1984) and the Geotiffs projected to the appropriate WGS UTM zone. However, when I loaded the indice and some of the photos together into ArcMap the images are significantly displaced to the north (whether I had the dataframe as geographic or projected). The indice is located in the correct position. Experimentally, I used ArcCatalog to clear all the transformation information. This resulted in the images lining up correctly with the indice, with their units still in meters. I didn't fundamentally like the idea of leaving the images without projection information, but:

Setting the projection to GCS WGS 1984 resulted in the images
disappearing completely. If I right click and choose "Zoom to Layer",
it leaves the window wherever it was at. The map extent also does no
change. The layer extent becomes "1,#QNAN", a null value.
Setting the projection to the WGS 1984 UTM results in the
original displacement (obviously)
Setting the projection to GCS WGS 1984 and projecting into UTM results in an unknown error (a 999999 error), which makes sense as there is an error in the aforementioned GCS projection.

I know this has probably  been asked a million times here in this forum and elsewhere on the web, but I am stumped.


texts: change detction for landsat-8 and landsat-5 images

texts:
Can anyone help me in suggesting methods for change detection for landsat-8 and landsat-5 images?

The software that I am using is ArcMap 10.2 and i need to find the increase in builtup area from the two images.

I tried it in ArcMap 10.2 using image analysis tools's difference tool but unfortunately it is not enabled.

My work is to get the increase in builtup area from landsat 2005 image and landsat8 2015 image.


texts: England - Non-metropolitan counties and the ONSPD (Postcode data)

texts:
I have a large amount of data records which have a postcode on them that I'm in the process of trying to visualise on an England map.

The first map view I'm doing is naturally regional view. This is fairly trivial because the ONSPD (Ordnance Survey) postcode data has a Government regional field therefore I can lookup where each postcode belongs.

The issue I'm having is doing the same for the non-metropolitan counties that belong to each region. It seems that the ONSPD data does not seem to have a direct mapping to them. It maps to counties but those counties do not cover the whole of England.

The shapefile I have for the counties looks as such :

As you can see, this fills the entire England map and as those of you who have experience with the English county system there are shapefiles for counties that do not cover the whole of England etc.

I tried to contact the Office for National Statistics but they were unhelpful thus far. They simply asked me to have a look at their resources. I explained how I had gone through their website and the different type of county map files. I've also looked at the documentation on the ONSPD but unable to find a county level field which every postcode in England is mapped to. The county mapping they have of the postcodes looks like this :

Not including the Greater London boundaries (it's boroughs) it does not map onto the above linked shapefile.

Hopefully my question has made sense and someone else who has tried to visualise data on a county level for England has had the same issue.


texts: Geoserver: Problems with external graphic URLs

texts:
In Geoserver, I try to save a SLD file containing this Rule:

<sld:Rule>
  <sld:PointSymbolizer>
    <sld:Graphic>
      <sld:ExternalGraphic>
        <sld:OnlineResource xmlns:xlink="http://www.w3.org/1999/xlink" xlink:type="simple" xlink:href="https://my.url.net/services/pictureservice/wmg.pic?reference=jvb&id=1047068"/>
        <sld:Format>image/png</sld:Format>
      </sld:ExternalGraphic>
      <sld:Opacity>0.9</sld:Opacity>
      <sld:Size>14</sld:Size>
    </sld:Graphic>
  </sld:PointSymbolizer>
</sld:Rule>

The vaildator complains like this:

  line -1: The reference to entity "id" must end with the ';' delimiter.

The problem seems to be the & in the url, which leads to the validation error. But unfortunatly, in the url i need both parameters, the reference and the id. Is there a way to get around this?
I am on Geoserver 2.7.2


texts: How to Select by Attributes on features sharing common value?

texts:
I am attempting to 'select by attributes' a point feature with the lowest numeric value which shares a common value with other point features in another field. For example: I may have 3 point features with address values: 123, 125 and 127. Each of these point features shares a parcel address value of 131 in another field.

I simply need to select the lowest address point for each parcel address so that I can use a definition query to display a subset of the data (Only points with these lowest values). Thank you in advance for the help.


texts: QGIS: Unable to add a new feature when joining two PostGIS tables

texts:
I have two Postgresql tables: parent and child, as the following:

CREATE SEQUENCE public.parent_id_seq;

CREATE TABLE public.parent
(
    id integer NOT NULL DEFAULT nextval('public.parent_id_seq'::regclass),
    type character varying COLLATE pg_catalog."default",
    CONSTRAINT parent_pkey PRIMARY KEY (id)
);

CREATE SEQUENCE public.child_id_seq;

CREATE TABLE public.child
(
    id integer NOT NULL DEFAULT nextval('public.child_id_seq'::regclass),
    name character varying COLLATE pg_catalog."default",
    parent_id integer,
    CONSTRAINT child_pkey PRIMARY KEY (id),
    CONSTRAINT child_parent_id_fkey FOREIGN KEY (parent_id)
        REFERENCES public.parent (id) MATCH SIMPLE
        ON UPDATE NO ACTION
        ON DELETE NO ACTION
);

In QGIS I have made a join between these two tables (as shown in the attached image). But I am not able to add new features from the attribute table. When I do so I got the following error:

Could not commit changes to layer child

Errors: ERROR: 1 feature(s) not added.

  Provider errors:
      PostGIS error while adding features: ERROR:  invalid input syntax for integer: "nextval('public.parent_id_seq'::regclass)"
    LINE 1: ...ild"("id","name","parent_id") VALUES ($1,'name 9','nextval('...

How to Reproduce

1- Create these two tables, as described above.

2- Add them to QGIS.

3- Right click on the parent table, then select properties, then from the Joins tab create a new Join as shown in the attached screen shot.

4- Then allow editing both layers.

5- After that open the attribute table of the parent layer, then select Add Feature, and fill in the type and name fields.

6- Save the parent layer, it will be saved with no problems. Then save the child layer, it will generate the error mentioned above.

I am using QGIS 3.10


texts: Export User Statistics from ArcGIS Online

texts: Is there any way to export usage statistics to csv from an ArcGIS Online institutional account?  I would specifically like to be able to identify the number of times logging in or total time logged in as well as number and type of items created.


texts:
I'm trying to find a point along a series of line at a specified percentage of each line (the percentage along the line changes with every point too). I found ArcGIS documentation of the data management tool GeneratePointsAlongLines, but when I try to use it I get the error message:

Traceback (most recent call last):
  File "L:\gathr\indonesia\Sara\Scripts\1234.py", line 11, in
    arcpy.GeneratePointsAlongLines_management(Input_Features, Output_features, "PERCENTAGE",percentage)
AttributeError: 'module' object has no attribute 'GeneratePointsAlongLines_management'

What am I doing wrong here? Is there an easier way to go about this rather than Generate Points Along Lines. I'm having a hard time figuring out how to use positionAlongLine, but it looks like it could be helpful. Thoughts?

import arcpy
import os
in_data = "L:\\gathr\\indonesia\\Sara\\Lines_withinTimes\\Day_andmmsi\\s20160101fc_lines.gdb\\Detection1_lines_UTM"
rows = arcpy.SearchCursor("L:\\gathr\\indonesia\\Sara\\Lines_withinTimes\\Day_andmmsi\\s20160101fc_lines.gdb\\Detection1_lines_UTM")
for row in rows:
    row_name = str(row)
    Input_Features = in_data
    Output_path = "L:\\gathr\\indonesia\\Sara\\Lines_withinTimes\\Day_andmmsi\\s20160101fc_lines.gdb\\"
    Ouput_features = os.path.join(Output_path, "Detection_" + row_name)
    percentage = row.getValue("duration_fraction")
    arcpy.GeneratePointsAlongLines_management(Input_Features, Output_features, "PERCENTAGE",percentage)


texts: QGIS export vector as CSV with Decimal Degrees (not Meters)

texts:
I am trying to convert a Satellite JP2 file to a CSV database showing the lat/lon and a column with a calculated NDVI value. Both my project CRS and layer CRS are set to EPSG::28992. I have done the following:

Add raster layer with the JP2
Use the raster calculator to calculate NDVI index
Use the Raster value's to points tool in the Processing Toolbox

When step 3 is completed the Identify Results panel is no longer showing the coordinates in Decimal Degrees but in meters (except for the values under '(Derived)'. The current map coordinates are still shown in Decimal Degrees. When saving the shapefile to CSV it's only possible to choose XY in meters.

Is there a way to export lat/lon coordinates?


texts: Automate Select and Export to PDF using Arcpy Search Cursor

texts:
I'm trying to automate a process to select each feature of a feature class one by one, zoom to it and export as pdf.  using the code below, the first feature is selected and the pdf is exported, and the loop stops there. How can I modify this code to loop through all records?

import arcpy

mxd = arcpy.mapping.MapDocument('CURRENT')
df = arcpy.mapping.ListDataFrames(mxd, "Layers") [0]
fc = "This_FC"

with arcpy.da.SearchCursor(fc, ("SHAPE@", "FeatureID", "OID@")) as cursor:
    for row in cursor:
        query_str = """ "FeatureID" = '{0}' """.format(row[1])
        arcpy.SelectLayerByAttribute_management(fc, "NEW_SELECTION", query_str)
        df.extent = row[0].extent
        df.scale = df.scale * 5

        arcpy.RefreshActiveView()
        arcpy.mapping.ExportToPDF(mxd, "C:/test/" + fc + "_" + str(row[1]) + "_" +  str(row[2]) + ".pdf")


texts: How to normalize elevation of GeoTIFF file data stored as unsigned integers?

texts:
I have a GeoTIFF file where the pixels are stored as unsigned 16 bit integers. I am using C# Gdal bindings to read the GeoTIFF file, using the Band.ReadRaster method, in order to produce an array of normalized values between 0 and 1.

I use ComputeRasterMinMax to compute the min and max value of all the data points in the band. For other GeoTiff files where pixels are stored as floats/doubles, it works out that the min and max are the absolute min and max elevation in meters. However, it is not clear to me if this is true when the pixels are stored as unsigned integers.

For instance, a sample file gives me 4000 as the min and 27000 as the max. If these are absolute elevations, then would I need to normalize each pixel via the method:

normalized pixel value = (pixel value - 4000) / 23000.

Is that correct, or should I normalize the pixels using ushort.MinValue and ushort.MaxValue, i.e.:

normalized pixel value = (pixel value - ushort.MinValue) / (ushort.MaxValue - ushort.MinValue)

which of course can be simplified given ushort.MinValue = 0.


texts: Connecting offset point symbol to its original position in QGIS?

texts:
I'm drawing a track where I want to place icons (point symbols), but with an offset so that the line is fully visible. Here is an example:
The origin is exactly on the track, therefore when zooming in and out the icons stay in a reasonable position. So far,so good.
What I need, though, is to connect the edge of the icon to the original point on the track using a line or arrow. I need this to indicate where exactly is the danger point (that is what the exclamation mark icon stands for)


texts: Using ArcGIS Tokens to visualize user authenticated data in Esri Leaflet?

texts:
I'm currently using Esri Leaflet with authenticated user services. Authenticated services are set up and work fine for displaying geographic data after login, however, I want to  create a dashboard and visualizations (mostly graphs) from the data in the SDE within our mapping platform. eg."Your trends over time vs all users trends over time". I can export the data as csvs from our database and use D3, however I'm really stuck on a more efficient way to do this using authenticated layers.

Is there a way to create visualizations while using an esri token to get the user-specific data? Additionally, I want to aggregate all data to display as well (from other users) but I am concerned about security. Eg. A graph of User A's submissions vs all submissions by  all users (user b, user c) without user b and user c's data being exposed to user A.

I can't find any working examples of this. I'm mostly concerned about security I do not want our users to be able to see other user's data (unless aggregated).

I've seen one similar unanswered question at Seeking data visualization techniques for web maps using Esri and open source suites?. This question addresses the D3 aspects but not authentication or querying the data.


texts: Create a 5x5 km polygon grid by sub-setting a 10x10km one

texts:
I have a 10x10km grid in Lambert Conformal Conical projection (consider it as a sample frame).

I want to create another grid 5x5km by splitting in 4 parts each one of the original grid cells (the idea is to decrease the sample unit size maintaining a certain degree of spatial comparability).

In ArcGIS I tried creating a 5x5km polygon fishnet using as template the original grid. However, probably due to the conical projection, the output has got a variable shift from the input.

I then found a "Custom Grid Tool" (please see for reference http://ianbroad.com/creating-quarter-quarter-section-grid-python/)
which creates a polyline to be used to split the original grid into its quarter section. Unfortunately the original grid has too many polygons and the splitting operation makes the pc to crash.

How can we do this using ArcGIS Desktop?


texts: How to clip a raster by extent using Processing from the Python console?

texts:
I want to clip a raster by extent in QGIS with the python console. Here is the description given in the documentation.
You have to use the following input parameter:

Input layer
Nodata value (no_data)
Clipping extent (projwin)
Additional creation parameter (extra)
output

The syntax is:

processing.runalg('gdalogr:cliprsterbyextent', input, no_data, projwin, extra, output)

I use the following data:

input = "C:/Users/....../inputraster.tif"

no_data= , , [I leave it empty]

projwin = 4278848, 3020676, 4297792, 3024646 [without " "]

extra = , , [I leave it empty]

output = "C:/Users/....../outputraster.asc"

I used these given examples in several variations, but I alway get a SyntaxError: invalid syntax.


texts: how to load in_memory data as input for an ARCGIS 10.1 model?

texts:
I'm trying to run a model that gets  stops from an in_memory workspace and creates a route with directions based on the route, here is what the model looks like:

The following error (error 030024: Solve returned a failure) Occurs:

I understand that the model doesn't recognize the stops I put in the in_memory, but as I couldn't find documentation on how to enter data into the in_memory and use it as model input (all documentation I found was how to use in_memory for the output), could anyone explain how it's done?

I need to make the stops come from in_memory (and not a shp) as the model is supposed to turn into  a geoprocessing package (gpk) in which the user will have to input the stops by clicking them on the map.


texts: Storing data from ArcGIS layer in SQL Server?

texts:
I'm  new to ArcGIS and I want to know if I can store data from a layer to a SQL Server table and if we can do this, then how?

I have already made a connection between ArcGIS and SQL Server.


texts: Counting the number on elements in a relation, with QGIS 3.x

texts: Given a relation in QGIS 3.x (I have a layer A of points (150 geometries) and an attribute table B of users with 500 records, and each one of this registered users is associated to a point geometry through a relation), is it possible to store in a new field of the layer A (points) the number of registered users that are associated to each point, using the established relation?


texts:
How can I use nested styles?

I found Nathan W.'s post "Improvements to the QGIS rule based rendering" about nested styles, but I do not understand how to use?


texts: Oracle add layer and Oracle Connection icons are missing

texts:
I have recently installed QGIS 2.8 and Oracle 11.2 Client on Linux Red Hat 7.

Everything was fine, I am able to access Oracle database, as well as running QGIS.
Unlike QGIS 2.8 from WIndows, Linux version does not have Add Oracle spatial layers and Oracle new connection icons.

I think the Oracle support was removed by QGIS. Is this correct?

The machine where QGIS runs has installed Oracle Client 11g, unixODBC and SQL Developer.

I have tested the oracle connectivity from a terminal windows as well as from SQL Developer and it works.

I tried as alternative to load spatial data the followings:

cx_Oracle module from python: it connects to oracle from QGIS python console but cx_Oracle does not have functions to add SDO date to a QGIS layer.
Add vector layer from Database ODBC- does not connect to an Oracle


texts:
I have to define a geometry column in PostGIS for 2D points, but when I try to import points from DXF to postgres I get the following message:

  "postgis error: when adding features geometry has z dimension but
  column does not"

Is there any way to change DXF files to be 2D geometry? I have tried this with the comand Flatten in autocad but the same error appears.


texts: Problem with massive data importation to PostGIS/PostgreSQL

texts:
I need to import a huge SIG data into PostgreSQL (my shapefile table have almost 1.700.000 rowns and 43 columns).
I'm using PostgreSQL 12.4 and PostGIS 3.0.x, that came together with PostGIS Bundle 3 DBF Loader, which I'm using for the importation.
The thing is that my data is being loaded, at first look with no problem, but verifying the integrity of the table, I noticed the information are all messed up - with column values misplaced between them (the number of rows is aparently ok).
Does someone know the reason for that? Is that a size limit for data to be sucessfully loaded from external environment into PostgreSQL/PostGIS?
I'm new in this topic, and don't know what other information is important for someone to help me.


texts: How to use geoserver to display data on hosted website

texts:
I am using geoserver to display GIS data on my hosted website
when geoserver is running I can see my data on my website in both firefox and Internet explorer but only on my laptop where I have geoserver running, when I try to view the website using Firefox or IE on another PC I just get a blank box with error on page warning on the bottom of the web browser. can anyone help with this thanks


texts: Convert a set of points to multiple single points?

texts:
I have a set of 560 building points and a set of many thousand points representing each meter of a tunnel that is under construction. I don't (only) need to know the nearest, but the distance between each building (red) and every single tunnel point (in blue):

Because I only have the Basic licence (with 3D analyst + Spatial analyst extensions) I plan to run the Near 3D tool with Input feature = Tunnel points and Near feature = One single building point. The problem is that I have 560 building points and I need to separate or split the point feature class into 560 single feature classes, each representing one building.

Is there an easy way of splitting a feature class with many points into many feature classes containing only one point?

I have started with selecting one point and then exporting it into a new feature class manually, but this will be a lengthy process (and I'll probably need to do it with more than the first 560 points).

I am not well-versed with python or scripting, but am open to such a solution.


texts: Any Good Map Rendering Engines for SqlGeometry/SqlGeography

texts:
Can anyone suggest good map rendering engines for MS SQL Server Spatial types (SqlGeometry/SqlGeography)?  By rendering, I mean plot the geometries to a common image format such as png or jpeg with very simple symbology.

The 2 I have experimented with thus far are:
1.  SharpMap
2.  MapServer with the MSSql plugin

I just want to see if there is anything else out there that I might have missed.


texts: change the area in a osm-request: from lon & lat to whole country

texts:
I am new to PHP's SimpleXML. the script below makes a starting point for my project. Below I have added my own answer to the question, how to refine the code to ad more tags.

I want to filter the data to get the nodes with special category. Here is sample of the OSM data I want to get the whole schools within an certain area.

The first script below runs well - but now I want to refine the search and add more tags. I want to store all into MySQL.So we need to make some XML parsing with PHP:

The following is a little OSM Overpass API example with PHP SimpleXML

NOTE: i want to change the values in the boundary-box. So that we can gather the shools from a certain area - eg. the country Italy

Q: is this possible?

see the script:

<?php
/**
 * OSM Overpass API with PHP SimpleXML / XPath
 *
 * PHP Version: 5.4 - Can be back-ported to 5.3 by using 5.3 Array-Syntax (not PHP 5.4's square brackets)
 */

//
// 1.) Query an OSM Overpass API Endpoint
//

$query = 'node
  ["addr:postcode"~"RM12"]
  (51.5557914,0.2118915,51.5673083,0.2369398);
  node
  (around:100000)
  ["amenity"~"school"];
           out;';

$context = stream_context_create(['http' => [
    'method'  => 'POST',
    'header' => ['Content-Type: application/x-www-form-urlencoded'],
    'content' => 'data=' . urlencode($query),
]]);

# please do not stress this service, this example is for demonstration purposes only.
$endpoint = 'http://overpass-api.de/api/interpreter';
libxml_set_streams_context($context);
$start = microtime(true);

$result = simplexml_load_file($endpoint);
printf("Query returned %2\$d node(s) and took %1\$.5f seconds.\n\n", microtime(true) - $start, count($result->node));

//
// 2.) Work with the XML Result
//

# get all school nodes with xpath
$xpath = '//node[tag[@k = "amenity" and @v = "school"]]';
$schools = $result->xpath($xpath);
printf("%d School(s) found:\n", count($schools));
foreach ($schools as $index => $school)
{
    # Get the name of the school (if any), again with xpath
    list($name) = $school->xpath('tag[@k = "name"]/@v') + ['(unnamed)'];
    printf("#%02d: ID:%' -10s  [%s,%s]  %s\n", $index, $school['id'], $school['lat'], $school['lon'], $name);
}

?>

if we change this

$query = 'node
  ["addr:postcode"~"RM12"]
  (51.5557914,0.2118915,51.5673083,0.2369398);
  node
  (around:100000)
  ["amenity"~"school"];
           out;';

$context = stream_context_create(['http' => [
    'method'  => 'POST',
    'header' => ['Content-Type: application/x-www-form-urlencoded'],
    'content' => 'data=' . urlencode($query),
]]);

we get more than 2170 results:

#2171: ID:2837001831  [51.1783749,-0.7970866]  More House School
#2172: ID:2865091022  [50.9090614,-0.5565425]  Dorset House School
#2173: ID:2882477853  [51.6261198,-0.7349665]  Bowerdean Primary School
#2174: ID:2901434856  [51.6542477,-0.3098923]  The Haberdashers' Aske's Girls School
#2175: ID:2901434857  [51.6565707,-0.3129822]  The Haberdashers' Aske's Boys School

NOTE: i want to change the values in the boundary-box. So that we can gather the shools from a certain area - eg. the country Italy or something thelike

is this possible?

update: can i do some rework of the area-data like so:

[place] => SimpleXMLElement Object
        (
            [@attributes] => Array
                (
                    [place_id] => 331526
                    [osm_type] => node
                    [osm_id] => 96209423
                    [place_rank] => 15
                    [boundingbox] => 44.4361381530762,44.4361419677734,26.1027431488037,26.1027450561523
                    [lat] => 44.436139
                    [lon] => 26.1027436
                    [display_name] => BucureČti, Sector 2, Bucuresti, RomĂ˘nia, European Union
                    [class] => place
                    [type] => city
                    [importance] => 0.73231672860554
                    [icon] => http://nominatim.openstreetmap.org/images/mapicons/poi_place_city.p.20.png
                )

        )

is this doable!`?


texts: Using split command' in Field Calculator?

texts:
I have bunch of home adresses which contain home number and subnumber, for example Street of J.F. Kennedy 1005A. Now I want to split this into two fields out of home number, one for main number and one for subnumber.

How to split '1005A' and get '1005' in one and 'A' in another field somehow using the 'split command' in Field Calculator?


texts: Instructions for building a Flex App using open source PostGIS DB and Geoserver

texts:
Anyone have a link for a good guide on how to build a Flex app to work with Open source REST services? I am familiar with how to do it in ESRI and Arc Server Manager, but don't want to pay for a ArcServer license.

Thanks!
-Adam


texts: Age of RTK corrections - fixed solution accuracy

texts:
Does anyone know how the age of NTRIP corrections influences accuracy of an RTK solution?

I read somewhere that correction loss for few seconds or even a minute may still be acceptable for some applications. I ask this because I live in Europe where we have a nice dense CORS network equivalent but most countries charge a fee for its usage based on time. So one could save quite some money by asking for corrections less frequently (maybe every 10 seconds or even less, given that payments units are seconds).

I found one post with vertical test (unknown agricultural-grade setup): http://lefebure.com/articles/rtk-correction-data-age-accuracy/

I am interested in horizontal accuracy for a L1 and L1/L2 GNSS RTK setup.


texts: Calculating distances between buildings, plots and road network

texts:
I have three layers, see image below:

A polygon layer containing building outlines (Purple)
A polygon layer containing plot outlines (Yellow)
A polyline layer showing the road network ( Red)

What I want to do is calculate whether there is room for a driveway on the plot.
Therefore I want to calculate the distance between every building and the plot edge (plot polygon containing all or part of a building)  that is closest to the road network.


texts: What is difference between grid layer and feature Layer in mapbox map?

texts:
I have created polygon using two way but i have not understand difference between them.

using grid layer

http://bl.ocks.org/rclark/7043524
using feature layer......

L.mapbox.accessToken = 'pk.eyJ1IjoiZ2F1cmF2Y2F0c3RlY2giLCJhIjoiY2l1cGo0MTl4MDAxajJ1bng5a2xieTY0diJ9.0fDXNulI91U85ngSc4jGCg';

var map = L.mapbox.map('map')
.setView([38.89399, -77.03659], 17);

L.mapbox.styleLayer('mapbox://styles/gauravcatstech/civ82vos7008f2is5x92cpvye').addTo(map);

var featureGroup = L.featureGroup().addTo(map);
var circle_options = {
    color: '#fff',      // Stroke color
    opacity: 1,         // Stroke opacity
    weight: 10,         // Stroke weight
    fillColor: '#000',  // Fill color
    fillOpacity: 0.6    // Fill opacity
};

var line_points = [
    [38.893596444352134, -77.0381498336792],
    [38.89337933372204, -77.03792452812195],
    [38.89316222242831, -77.03761339187622],
    [38.893028615148424, -77.03731298446655],
    [38.892920059048464, -77.03691601753235],
    [38.892903358095296, -77.03637957572937],
    [38.89301191422077, -77.03592896461487],
    [38.89316222242831, -77.03549981117249],
    [38.89340438498248, -77.03514575958252],
    [38.893596444352134, -77.0349633693695]
];

var polyline_options = {
    color: '#000'
};

var polyline = L.polygon(line_points, polyline_options).addTo(featureGroup);

var drawControl = new L.Control.Draw({
    edit: {
        featureGroup: featureGroup
    }
}).addTo(map);

map.on('draw:created', function(e) {
    featureGroup.addLayer(e.layer);
});


texts: Adding non geographic data from MSSQL Database to QGIS?

texts:
In QGIS (2.8) I am trying to add non-geometric data / flat data that originates in a MS SQL database. This flat data later should be joined with other layers etc.

Connecting to the specific database = working
Listing the table that I would like to open up in QGIS = working

When adding the selected 'flat' table I get the following error message:

dbname='data' host=server\sql2008 table="dbo"."TABLE" sql= is an invalid layer - not loaded

Is there a "trick" or a set of conditions that the table needs to meet before it can be added into QGIS as "flat" data?


texts: "Couldn't load PyQt" While starting QGIS 3.0

texts:
The following error shows when I am trying to open QGIS 3.0. Although I have installed PyQt5 and the module can be imported from python.

I am using Ubuntu 16.04 LTS


texts: How to render a .tiff file on Gmaps by using Cloud Optimized GeoTIFF server?

texts:
I am doing a Web project where I have used GMaps to render a (.tiff) file which I converted to Google Map Tiles. Now I want to do this dynamically where the user just uploads the (.tiff) file and it shows on Gmaps. I have come across COG (Cloud Optimized GeoTIFF) server and I see that might be able to solve my problem. But I am not getting to understand the general workflow of the whole process and to set up a COG with support for GMaps as a base map.
How does COG work?


texts: Using IF/ELSE IF/ELSE statements in Attribute Expressions for String Fields

texts:
I am trying to do build some Attribute Expressions using Arcade and I am having some issues with it returns the else statement. I know that there are features that would return true because I have used these parameters under filters before.

var parksign = $feature.ParkSign2;
var signpole = $feature.SignPole1;

if(parksign == 'One on bus stop pole, 2nd on breakaway in ground, grass, dirt'){
    return 'Yes'
}else if(signpole == 'Yes, Metro breakaway in ground, grass, or dirt'){
    return 'Yes'
}else{
    return 'No'
}


texts: Setting layer order for output in QGIS graphical modeler?

texts:
I've recently discovered the graphical modeler in QGIS and I've managed to speed up my workflow significantly when doing ecological and geo/hydrological analyses on a landscape level. Al though I still haven't figured out how to get my output layers to stack as I want them to.
A typical model looks like this (some Swedish in there):
I input a DEM and two kinds of landscape erosion maps (one vector and one raster) and the model cuts them from a nation wide extent down to the extent of the DEM and then styles them. The model also creates a hill-shade layer and contours.
The "problem" arise when the model completes and and all my newly created layers end up scrambled and I have to order them by hand (yes, I'm lazy) - is there a way to get them to come out in the order I want: hill-shade at the bottom with the erosion maps on top and finally the contours?


texts:
I am using the plugin Leaflet DVF (https://github.com/humangeo/leaflet-dvf). It uses this code to bind a popup to a polygon and create a table inside it displaying its corresponding records (a row for each record) coming from a JS file.

This is the code for the popup:

                onEachRecord: function (layer,record) {
        var $html = $(L.HTMLUtils.buildTable(record));

        layer.bindPopup($html.wrap('
').parent().html(),{
            minWidth: 600,
            maxWidth: 700
        });
    },

This is the file with the record information:

var NLpopdens =[
{
    "CODIGO": "01",
    "NAME": "10",
    "POPU": "635100",
    "AREA": "162.18",
    "Density": "Title
"
},
{
    "CODIGO": "02",
    "NAME": "11",
    "POPU": "263500",
    "AREA": "227.11",
    "Density": "Title

"
},

I am able to run the popup with an image but not with a DIV inside it. Does anyone know how to solve this?

PS. The Console in the browser gives me the following error:

  Uncaught Error: Highcharts error #13: www.highcharts.com/errors/13
      at Object.a.error (highcharts.js:10)
      at a.Chart.getContainer (highcharts.js:258)
      at a.Chart.firstRender (highcharts.js:273)
      at a.Chart.init (highcharts.js:249)
      at a.Chart.getArgs (highcharts.js:248)
      at new a.Chart (highcharts.js:247)
      at Object.a.chart (highcharts.js:247)
      at HTMLDocument. (conflictdata.html:679)
      at c (jquery-1.9.1.min.js:3)
      at Object.fireWith [as resolveWith] (jquery-1.9.1.min.js:3)

And the Highcharts code:

Highcharts.chart('container', {
chart: {
    renderTo: 'container',
    type: 'area'
},
title: {
    text: 'Area chart with negative values'
},
xAxis: {
    categories: ['Apples', 'Oranges', 'Pears', 'Grapes', 'Bananas']
},
credits: {
    enabled: false
},
series: [{
    name: 'John',
    data: [5, 3, 4, 7, 2]
}, {
    name: 'Jane',
    data: [2, -2, -3, 2, 1]
}, {
    name: 'Joe',
    data: [3, 4, 4, -2, 5]
}]

});


texts:
I want to use arcpy raster calculation to create some script but the output raster to get from this arcpy raster calculation script is different from this to take from model builder in ArcGIS gui,any idea why ?

arcpy script :

con='("{0}"*"{1}"*"{2}")'.format(raster1,raster2,raster3)
arcpy.gp.RasterCalculator_sa(con, rstercalc)

raster calculation in model builder :

my raster is float if I try to use float(raster1) I take error :

float cannot convert in memory layer

I have some wrong in con in script or the wrong is the in_memory layer where I use ?


texts: Performing Viewshed Analysis in ArcGIS Desktop

texts:
I'm wondering if anyone can help me out. I'm trying to carry out a viewshed analysis of wind farms and single turbines on a region of the UK. I have the landform data as well as the coordinates of around 1,600 turbines, I'm looking to create a 3D map, plot these turbines on the map and carry out a viewshed analysis to determine what area is suffering from the highest level of visual impact.

I am using ArcGIS Desktop 10.2.1.


texts:
I have a leaflet map that consists of various wms layers representing different years.

I am trying to add a legend to my leaflet map but with my own custom hex colours instead of trying to use a layer in the map (as my legend is more broad) This is the sample code from leaflets website.

var legend = L.control({position: 'bottomright'});

legend.onAdd = function (map) {

    var div = L.DomUtil.create('div', 'info legend'),
    grades = [0, 10, 20, 50, 100, 200, 500, 1000],
    labels = [];

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

return div;
};

legend.addTo(map);

Is there anyway to adapt this to my own 5 class hex(colors) legend?


texts:
I used ArcGIS javascript API 3.0 version and I use Bing map. I need to be able to zoom to 22., 23., 24., .. levels. So I was added an addNumLevels function to Map;

 esri.Map.prototype.addNumLevels =function(num){

    for(var i = 0; i < num; i++ ){
        var ln = this.__tileInfo.lods.length;
        this.__tileInfo.lods[ln] = {};
        this.__tileInfo.lods[ln].level = ln;
        this.__tileInfo.lods[ln].resolution = this.__tileInfo.lods[ln-1].resolution / 2;
        this.__tileInfo.lods[ln].startTileRow = this.__tileInfo.lods[ln-1].startTileRow;
        this.__tileInfo.lods[ln].startTileCol = this.__tileInfo.lods[ln-1].startTileCol;
        this.__tileInfo.lods[ln].endTileRow = this.__tileInfo.lods[ln-1].endTileRow * 2;
        this.__tileInfo.lods[ln].endTileCol = this.__tileInfo.lods[ln-1].endTileCol * 2;
    }
}

But I when I upgraded the ArcGIS js api to version 3.2, my addNumLevels function stopped working.
Can anybody help me?

Thanks.


texts:
I am trying to clip a raster image with the outline shape of a country, I imported the .png image in using the georeferencer plugin. I then selected the country I wanted to use and saved it as a shapefile. I then tried to use the Raster extraction - clip raster by mask layer.

But this did not work

The last image is what I managed to get it to output, but on the other occasions it did not work at all and that last image is not exactly what I want


texts: ArcSDESQLExecute: Mismatching number of BEGIN and COMMIT statements

texts:
I have the following sql update query that I am using as part of a larger python script.   It essentially takes data from one Feature Class (FC) to  update a Versioned View of a different FC.

Important to note the SQL query runs with no errors in SQL Server Management Studio but it fails when trying to run in python.

Code:

#Set Variable for SDE connection
# ArcSDESQLExecute allows SQL statements to be run inside the Python script
sdeConn = arcpy.ArcSDESQLExecute(Environ)
sql_statement = '''EXEC devsde.sde.set_current_version "SYOKLIC1.test_1";
       EXEC devsde.sde.edit_version "SYOKLIC1.test_1", 1;
       BEGIN TRAN;

       UPDATE [devsde].[syoklic1].[BOOSTER_evw]
       SET [devsde].[syoklic1].[BOOSTER_evw].[MaintenanceArea] =
          (
            SELECT [devsde].[syoklic1].[BOOSTER_Join].[MaintenanceCode]
            FROM [devsde].[syoklic1].[BOOSTER_Join]
            WHERE[devsde].[syoklic1].[BOOSTER_Join].[VMID]
            = [devsde].[syoklic1].[BOOSTER_evw].[VMID]
            AND [devsde].[syoklic1].[BOOSTER_Join].[MaintenanceArea]
            != [devsde].[syoklic1].[BOOSTER_Join].[MaintenanceCode]
          );

         COMMIT;
         EXEC devsde.sde.edit_version "SYOKLIC1.test_1", 2;'''
sql_statement_list = sql_statement.split(";")
for query in sql_statement_list:
    try:
        # you first need to designate which transactioanl version sql should be inteacting with
        sdeConn.execute(query)
        print "execution started " + "{{{" + query + "}}}"
    except Exception as exc:
        print "execution failed " + "{{{" + query + "}}}"
        print traceback.format_exc()
        print exc
        pass

The error comes up for the BEGIN TRAN and COMMIT statements

Error text:

ArcSDESQLExecute: StreamExecute ArcSDE Extended error 266 [Microsoft][SQL
Server Native Client 11.0][SQL Server]Transaction count after EXECUTE
indicates a mismatching number of BEGIN and COMMIT statements. Previous
count = 1, current count = 0.

How do I work around this?

Perhaps a way to format so that ArcSDESQLexecute can read the entire statement where the number of Begins/Commits are equal.


texts: Mass downloading raster tiles using QGIS

texts:
I want to mass download georectified sat image tiles (XYZ tiles) and save them locally.
I could setup QGIS so I can see google sat images.
So for example, I want all tiles of Los Angeles in a geo resolution of 50x50 m (and let the pixel resolution be whatever it will be then) and 200 x 200 tiles of that.
How I am going to approach this problem?


texts:
When publishing a Web Map Image Layer from ArcGIS Pro where the map uses a font in the symbology, I get this warning:

  20025 Layer uses fonts which are not embeddable and are not on the server

This would normally be related to the font is not installed on the server I try to publish to, and the first time I tried to publish, the font in question was not installed on the server.

But now I have installed the font on the server but I still get this warning!

I ignored this and ran the publishing anyway. It went smoothly, and in the resulting published service, the symbols are shown correctly!

So why do I still get this warning in ArcGIS Pro?


texts: Coded value domains, Allow Null Values = no. Why does it default to zero?

texts:
Trying to use Coded Value Domains with Allow Null Values = No in a File GDB, ArcGIS 10.2. Through an GIS.SE question that I can't seem to find again, I found this through ArcGIS help:

"Note:
The geodatabase model is such that it will insert an empty value (numeric = 0, text = "") instead of a database NULL if, and only if, the field in the database has a NOT NULL constraint on it."

I get that this prevents Null values from being entered. But it defaults to blanks or zeros, which overrides my domains when a new feature is created! So someone can add a feature and accept 0 in a field that has a coded value domain, which to some degree defeats the purpose of having the Coded Value domain to force certain values.

What I would like is to force someone to select a value before it will allow the feature to be saved. I thought this was the point of a domain?

Am I missing something or not understanding something about databases? Is this just a limitation, or am I trying to get it to do something that it inherently shouldn't do?


texts: Changing 'QgsVectorLayer' Attribute (QuickWKT)

texts:
I'm modifying QuickWKT plugin by elpaso, which is available on github: https://github.com/elpaso/quickwkt/blob/master/QuickWKT.py

This plugin stores 'layer' data inside memory layer (scratch memory)

What I'd like to do is storing the 'layer' data for later use (Qgis project) instead of scratch memory, so user can input new coordinates and generate geometry through the plugin

When I changed the vector data providers atributes of QgsVectorLayer from 'memory' (Line:121):

layer = QgsVectorLayer(typeString, layerTitle, "memory")

to 'ogr':

layer = QgsVectorLayer(typeString, layerTitle, "ogr")

After changing to 'ogr', the program doesn't output anything when running QuickWKT (see screenshot). I'm not sure whether the probem was caused from QgsVectorLayer's attribute doesn't support ogr, or the code block which invokes 'layer' doesn't support ogr.

Is it possible to change the vector data providers atributes from 'memory' to something else?


texts:
I am trying to create a user-defined defined function for a LAScatalog object in R using the lidR package and no matter which function I create, if opt_wall_to_wall is set to TRUE, the code returns the following error:
Error: User's function does not return NULL for empty chunks. Please see the documentation of catalog_apply.

Here's how I am reading in my LAScatalog:
LidRtst <- lidR::readLAScatalog(folder = "", filter = "-keep_class 2 3 4 5 8 9 26 -drop_x_below 7648119.63 -drop_x_above 7653883.38 -drop_y_below 675298.33 -drop_y_above 681872.34 -drop_overlap -drop_synthetic -drop_withheld -drop_number_of_returns 0 -drop_z_below 0")
lidR::las_check(LidRtst)
lidR::projection(LidRtst) <- rgdal::showWKT(as.character(LidRtst@proj4string))
lidR::opt_wall_to_wall(LidRtst) <- TRUE
lidR::opt_chunk_buffer(LidRtst) <- 100
lidR::opt_progress(LidRtst) <- TRUE
lidR:::catalog_laxindex(LidRtst)
future::plan(strategy = future::multisession, workers = detectCores())
lidR::set_lidr_threads(detectCores())

And here's a test to see what happens if I create a user-defined function with wall-to-wall set to TRUE
test = function(cluster) {
  las = readLAS(cluster)

  if(is.empty(las)){
    return(NULL)
  } else {
    return(list(0))
  }
}

opt_output_files(LidRtst) <- paste0(tempdir(), "/file_{ID}")
ret = catalog_apply(LidRtst, test)
#> Error: User's function does not return NULL for empty chunks. Please see the documentation of catalog_apply.

Now if I set wall to wall to FALSE
lidR::opt_wall_to_wall(LidRtst) <- FALSE
ret = catalog_apply(LidRtst, test)
head(ret)

It works! So how do I get a user-defined function in lidR to work wall-to-wall?


texts: Formula to create representative-sized rectangular map cell grids from degrees

texts:
I'm a novice user in R.

I have a dataset consisted of coordinates of some "kind of centroids" of grid cells and my question is if any formula exists to create representative sized grid cells using the size parameter of geom_point, knowing the size of my cell in degrees.

Here is an example of my code:

x <- ggmap(map) +
      geom_point(data = df,
                aes(colour= zscore), shape=15,alpha=0.8,

                size = varsize)

For example, in case I know the cell grids are 0.05 in degrees, which size would be preferred to use in my code ?

Here is an example of my grid using size=2


texts: Store geometry change on modify interaction OpenLayers

texts:
I am using OpenLayers modify interaction with events modifystart and modifyend. I am trying to remember state of feature before it is modified.
I do it this way:

var modificirane_cestice=new Object()

modify.on('modifystart',function(e){
    e.features.forEach(function(feature){
        if(modificirane_cestice[ogc_fid]){
            return }
        else{
            modificirane_cestice[ogc_fid]=new Array()
            modificirane_cestice[ogc_fid].push(new ol.Feature({geometry:feature.getGeometry()}))
            }

    })
    })

 modify.on('modifyend',function(e){
    e.features.forEach(function(feature){
    izmjenjena_geometrija=feature.getGeometry()

    modificirane_cestice[ogc_fid].push(new ol.Feature({geometry:izmjenjena_geometrija}))
    posljednja_geometrija=modificirane_cestice[ogc_fid].length-1
    pg_pov=ol.sphere.getArea(modificirane_cestice[ogc_fid][0].getGeometry())
    zadnja_pov=ol.sphere.getArea(modificirane_cestice[ogc_fid][posljednja_geometrija].getGeometry())
    console.log(pg_pov+".... zadnja pov...."+zadnja_pov)
    postotak=pg_pov/zadnja_pov
    console.log(postotak)
    if(postotak<0.9){
        console.log("postotak prelazi 10% i iznosi: "+(1-postotak)*100+"%")
    }
    })
    })

I expected to get two different areas and features but in the end I got the same one. I know there is problem in referencing feature object but don't know how to solve it. I want to create a new object which are independent on their source object.


texts: Choosing Java web framework for project?

texts:
We are developing java web gis application.

We use geoserver, postgis, geotools(all the core code).

We need to show/hide(or create/delete) dynamically vector layers(from shape file) and raster layers(wms layer from geoserver).
Also, our shapefiles are colorized(there are about 5-10 colorized zones in one shapefile).

Typical usage: user check/uncheck checkboxes -> layers are showed/hid.

We used geomajas in our 0 version, but there is the critical bug for us( http://jira.geomajas.org/browse/RTC-2 ).

So, could you please suggest any other java based web frameworks that can cover our needs?


texts:
I would like to run an external executable written in fortran which provides one of the inputs inside a python code.

Is it possible to run that exe inside Python script?


texts: Checking if point is inside area using OpenLayers?

texts:
I have a map and one point. I want create an area and check if a point is inside the area.

map = new OpenLayers.Map("map");
map.addLayer(new OpenLayers.Layer.OSM());

var lonLat = new OpenLayers.LonLat( -0.1279688 ,51.5077286 )
      .transform(
        new OpenLayers.Projection("EPSG:4326"), // transform from WGS 1984
        map.getProjectionObject() // to Spherical Mercator Projection
      );

var zoom=16;

var markers = new OpenLayers.Layer.Markers( "Markers" );
map.addLayer(markers);

markers.addMarker(new OpenLayers.Marker(lonLat));

map.setCenter (lonLat, zoom);

var bound = new OpenLayers.Bounds(5, 42, 10, 45).toGeometry();


texts: R-ArcGIS Bridge - Interpolation with gstat library

texts:
I have a problem with the idw or krige function from the gstat library in R. I'm trying to use the gstat library and the R-ArcGIS bridge to load a shapefile from ArcGIS into R, interpolate the data and then transform it back again. In the end I want to get a TIFF file with the interpolated grid, which I can load into ArcGIS.
Here is my script so far:
library(arcgisbinding)
library(sp)
library(gstat)

arc.check_product()

###--GET EXTENT AND COORDINATES FROM SHAPE--###
elev_shp <- arc.open("filepath/elev.shp")
elev_shp@extent
x_min <- 316116
y_min <- 5657570
x_max <- 317766
y_max <- 5659070

###--INPUT: LOAD DATA--###
elev_shp <- arc.select(arc.open("filepath/elev.shp"))
class(elev_shp)

###--CREATE GRID/RASTER--###
grd <- expand.grid(RW=seq(from=x_min, to=x_max, by=10), HW=seq(from=y_min, to=y_max, by=10))

###--CONVERT GRID TO SPATIAL OBJECT--###
coordinates(grd) <- ~RW+HW

###--ADD REFERENCESYSTEM TO GRID--###
proj4string(grd) <- "+proj=utm +zone=32N +datum=WGS84"
gridded(grd) <- TRUE
class (grd)
plot(grd)

##--CONVERT GIS-DATA TO SP_DATA--###
elev_shp_df <- arc.data2sp(elev_shp)
class(elev_shp_df)

bubble(elev_shp_df, zcol = "height", maxsize = 1, main = "Elevation by Points")

###--START INTERPOLATION--###
elev.idw = idw(elev_shp_df$height~1, elev_shp_df, grd)

##--CONVERT SP-DATA TO GIS-DATA--###
interpolation <- arc.sp2data(elev.idw)

###--OUTPUT: WRITE INTERPOLATION RASTER--###
arc.write("filepath", interpolation)

But when executing the function idw() I get the following error:
> elev.idw = idw(elev_shp_df$height~1, elev_shp_df, grd)
[inverse distance weighted interpolation]
Warnmeldungen:
1: In proj4string(d$data) :
  CRS object has comment, which is lost in output
2: In proj4string(newdata) :
  CRS object has comment, which is lost in output
3: In proj4string(newdata) :
  CRS object has comment, which is lost in outpu

Does someone have an idea of how I can solve this problem?


texts: Selecting/creating polygons of interest based on existing polygon database and indicate other polygons that share boundaries in ArcGIS for Desktop?

texts:
I am using ArcGIS 10.2 for Desktop.

The existing polygon database contains attribute table with the following structure, and has three levels; country, state, and city.  I am interested in the state of 'Anhui' for example (and many other states not listed below), and I want to create a field indicating states that share boundaries with any of them.

ID_0    country ID_1    state   ID_2    city    ID_state    ID_city
48      China   01      Anhui   08      Hefei   4801        480108
48      China   01      Anhui   08      Hefei   4801        480108
48      China   01      Anhui   15      Tongling4801        480115
48      China   02      Jiangsu 14      Suzhou  4802        480214
48      China   03      Fujian  15      Fuzhou  4803        480315
.
.


texts:
I am converting an .osm file to a .shp file with OGR2OGR:

ogr2ogr.exe --config OSM_USE_CUSTOM_INDEXING NO -skipfailures -f "ESRI Shapefile" output.shp input.osm

In osmconf.ini the other_tags is set to "yes" (other_tags=yes).

This creates an attribute other_tags which consists of a huge block of text, containing a lot of different tags.

Is it possible to somehow tell the OGR2OGR to split the other_tags into all of those tags it contains?

I searched, and found out that both in ArcGIS and QGIS there are options on how this can be done. but I am interested in doing it from OGR2OGR only.


texts: Scaling postgresSQL database when searching polygons against points?

texts:
I am using our single server -

CentOS 6.6 x64
16 CPU
24 GB
SSD
(PostgreSQL) 9.4.1

as our GIS database processing engine. This is part of the PostGIS software package installed on the server.

I am currently using near realtime polygon data against U.S. Census information to compute statistics. It works well and performance is acceptable.

What I would like to do, is add the ability mark homes/businesses and when a future polygon intersects that location, notify someone.

I have a sense of how I might program the SQL to do this, but I am really concerned about the performance. I run this query every hour on the previous hour of polygons I have received.

They are small polygons typically, but if I accumulate thousands of point locations to monitor, my concern is how I deal with performance.

I saw that the latest release of postgresSQL 9.6, has some support for parallel processing. I don't know if that would help in anyway.  I have also briefly read about some ability for postgresSQL to work in a cluster relationship across multiple servers.

Does anyone have a recommended strategy for building this the "right way" so that I do not run into performance issues and can scale the application?


texts:
I followed this tutorial and have some troubles. When I restart and get status of the renderd service, I get the following error :

renderd[8162]: An error occurred while loading the map layer 'default': Postgis Plugin: missing <table> parameter  encountered during parsing of layer 'landcover-low-zoom' in Layer at line 755 of '/home/anne/src/openstreetmap-carto/style.xml

It doesn't seems anyone else got the same error as me, I double-checked my renderd config file and the path to mapnik plugin seems to be correct :

user@comp:~$ mapnik-config --input-plugins
/usr/local/lib/mapnik/input

And in the renderd file:

[mapnik]
plugins_dir=/usr/local/lib/mapnik/input
font_dir=/usr/share/fonts
font_dir_recurse=true


texts:
I have an old ArcGIS 10.1 in my laptop. While I'm trying to create a geodatabase it gets created, but doesn't let me create new featureclass and/or feature dataset. I get the error:" failed to connect to database. bad login user."

Is there something I'm missing here?
I'm the admin on my laptop, and SQLExpress(V.2012), which has been installed by me. Though I don't know if I'm logging in as SQLEXPRESS' admin. I don't remember the credentials, I have definitely left it as default though.


texts: CartoDB, add search box after cartodb.createLayer

texts:
I used to use a hack I found online to add a search box to a cartoDB map, which I created using cartodb.createLayer and it used to work, but now it stopped working, the search box does not appear (no error message). Is there another way to do it?

.on('done', function(layer) {
      //hack from https://gist.github.com/javisantana/4e50d53e4eeb7fbdab0e
      var v = cdb.vis.Overlay.create('search', map.viz, {})
      $('#map').append(v.render().el);
    })

Full code is here: https://github.com/childcaremap/NYCdaycare/blob/gh-pages/dev/scoremap.html

Thanks!


texts: import .mbtiles[raster+png/jpg] from device can't view in Maptiler App

texts:
I am working with maptiler application in android[Kotlin+MVVM]. everything working fine except loading local .mbtiles map file. when I download a .mbtiles file from internet and load in to my android phone and then try to import it to my maptiler application, it will load in to database successfully but can not seen on screen . means there is some conflict to view that .mbtiles [raster+png/jpg] file. it gives following error:
Unable to parse resourceUrl maptiler://{filepath}

when remove maptiler:// from scheme and use direct file Uri it gives another error like:
Failed to load source maptiler-overlay-0: invalid value .at offset 0

now i did some try like change style.json with [style-local.json, style-cdn.json, style-mb.json] but cant get any success.
can Anyone had face this type of error or any idea how to resolve it ??
OR
can Anyone tell me the proper way to load local .mbtiles[raster+png/jpg]file in to maptiler[mapbox] app .??


texts: Why does arcpy.ListFeatureClasses() take so incredibly long?

texts:
I'm having an issue where calling arcpy.ListFeatureClasses() in a workspace with 400 shapefiles takes a very long time (more than 10 minutes).

Using the os module to achieve the same effect takes less than a second:

import os
allSingleNonRoadPolys = []
for file in os.listdir(r'D:\data\allSingleNonRoadPolys'):
if file.endswith(".shp"):
    print(file)
    allSingleNonRoadPolys.append(file)

You can then achieve the same effect as ListFeatureClasses by simply prepending the workspace to the filename, each time you need to use a dataset, as seen here:

for shp in allSingleNonRoadPolys:
    arcpy.someTool(workspace + '\\' + shp)

Why would I even keep using ListFeatureClasses() and what makes it take so long? Does it do anything incredibly useful that cannot be achieved in any other way?


texts: How to find out to which polygon each points belongs to in ArcMap

texts:
I have got a shapefile consisting of a number of polygons, and another shapefile with a number of points.

I want to add a column in the points shapefile. For each point, I want the column to contain the name of the polygon that surrounds it.

In simple terms, I just want to find out to which polygon each points belongs to. Any tips?

Here is the screenshot:


texts: DLL problem with QGIS 3.0 64bit (while 32bit works)

texts:
I installed the new version of QGIS 3.0 and first of all it worked all fine. After a while I tried to launch QGIS again and now the following problem appears.

After clicking "OK" another window appears:

So I opened the .env-file and it says:

    PATH=C:\OSGEO4~1\apps\qgis\bin;C:\OSGEO4~1\apps\Python36;C:\OSGEO4~1\apps\Python36\Scripts;C:\OSGEO4~1\apps\qt5\bin;C:\OSGEO4~1\apps\Python27\Scripts;C:\OSGEO4~1\bin;C:\WINDOWS\system32;C:\WINDOWS;C:\WINDOWS\system32\WBem
GDAL_DATA=C:\OSGEO4~1\share\gdal
GDAL_DRIVER_PATH=C:\OSGEO4~1\bin\gdalplugins
GDAL_FILENAME_IS_UTF8=YES
JPEGMEM=1000000
OSGEO4W_ROOT=C:\OSGEO4~1
PROJ_LIB=C:\OSGEO4~1\share\proj
PYTHONHOME=C:\OSGEO4~1\apps\Python36
QGIS_PREFIX_PATH=C:/OSGEO4~1/apps/qgis
QT_PLUGIN_PATH=C:\OSGEO4~1\apps\qgis\qtplugins;C:\OSGEO4~1\apps\qt5\plugins
QT_RASTER_CLIP_LIMIT=4096
VSI_CACHE=TRUE
VSI_CACHE_SIZE=1000000

What I have done so far:
1. uninstalling the program and reinstalling several times via Standalone-installer aswell as the OSGEO4W-installation - also tried the three different server and with Administrator-rights (as described in Problem starting QGIS 2.14, Prozedureinsprungpunkt is not found)
2. uninstalled and cleaned the registry-data as described in  How to completely uninstall and remove QGIS 2.8.1?
3. I searched for a duplicate of the Qt5Sql.dll that appeared another time in my Teamspeak-folder > so I uninstalled Teamspeak
4. Searched for leftovers from old versions and deleted all
5. uninstalled the QGIS-version via the OSGEO4W-installer and deleted the OSGEO4W-folder
6. installed the 3.0 (32-bit version) and later the 2.18.18-LTR version (no problems appeared there)

If anyone has another idea how to solve the problem, I am happy for your help!


texts: Land Use layer by combining, removing attributes and assigning new attributes

texts:
I'm working in ModelBuilder and I have limited (let's face it, no) Python knowledge, so I'd like to stay here if possible!

I have numerous shapefiles that I want to combine into a single land use layer, but I want the shapefiles to be classified into one of three categories. I'm not interested in keeping the original attributes; I just want the areas which can be combined and don't have to be separate.

My original files look like this:

A.shp - 90 rows, 10 fields
B.shp - 2500 rows, 15 fields
C.shp - 150 rows, 10 fields
D.shp - 500 rows, 3 fields
E.shp - 750 rows, 3 fields
F.shp - 150 rows, 15 fields
G.shp - 200 rows, 15 fields

and I want the final output to look like:

FINAL.shp - 3 rows (land use type A, B and C), 1 field (land use type)

The land use type is not specified anywhere in the original tables; I would have to specify this somewhere in the model.

If I were doing this manually in ArcMap, I would do something like:

Go through the shapefiles one by one, merging the rows.
Open the attribute tables one by one and delete the columns.
Create a blank shapefile and copy the land use A areas then merge them, then do the same for B and C.

Is there an easier way?


texts: How to find nearest Lat/Lon to the given Lat/Lon from set of Lat/Lon values?

texts:
I have a set of lat/lon values and gps lat/lon value. I want to find nearest lat/lon value to the gps point from that lat/lon set.

Is there any some formula or some other way so that i could implement that in my c++ program.


texts: Creating union of polygons with multiplying attribute table value in QGIS Modeler

texts:
I am new to QGIS.

I want to do union of two polygons A and B in such a way that one of the field values in attribute tables of A and B should be multiplied in the resultant union polygon C.

How can I do this in QGIS Modeler in one go?


texts: Displaying raster data by clicking on a field in an vector attribute table using QGIS 3

texts:
I am doing the following using QGIS 3:

Georeferencing of raster files.
Digitalization of polygons (which are in this raster).

It is often necessary to check how well the raster was digitized. To do this, I have to find the raster on disk and add it to QGIS.
In order to reduce the time required to search for a raster on disk (and manually add it in QGIS), I would like to put the raster data in an attribute field and, by clicking on this field, turn the raster on and off.
How can I do that?


texts: Calculating GPS coordinates given distance between points

texts:
I have 2 points where I know the GPS coordinates of one of the points. Additionally, I know their coordinates in some reference frame where I can calculate the distance between them in real world units (like meters).
So, the setup could be as follows
Point A: (GPS coordinates), (x1, y1)
Point B: (x2, y2)

where the simple Euclidean distance between A and B gives the distance between them in meters. For this, we assume that these points are fairly localized (spatially very close) and they both lie on the same plane.
Given this setup, could we calculate the GPS coordinate of point B?


texts:
I'm trying to get the newest version of Postgis (2.2.0) up and running with Postgres 9.4.

I'm using Windows and the only installation manuels I found, are Linux-only.

Can anybody point me in the direction of a useful Windows-installer?


texts: Vector Tiles VS Raster Tiles for leaflet

texts:
I have downloaded many vector tiles in the mbtiles format from OpenMapTiles. The reason I did that is because I will be using them on a station that cannot have an internet connection at all.

I would like in the end to see these tiles with leaflet on a browser, but it seems that the plugin for leaflet only considers raster-based tiles and not vector-based tiles, or so I thought until I saw Loading vector tile layer in Leaflet map? here on GIS. However, it seems to be outdated and the JSFiddle doesn't work anymore... So I figured I would try it out myself, but it did not work. So I am currently inclined to think that this plugin only works with raster tiles and not vector ones.

I also considered then converting my tiles to PNG so that I could use them with leaflet and used mbutils for that, as suggested on GIS at Finding good Mbtiles to png tiles converter? but obviously none of the generated PNGs are readable by my OS or leaflet as mbutils doesn't raster them.

I have two questions depending on how you  advise me:

Is it possible to read vector tiles with leaflet?
If not, how is it possible to convert them to raster?

I cannot use commercial services like maptobox or any similar services, this the reason why I was trying to have everything running on my machine right away.

I have seen before the vector grid from leaflet. It says that it can read vector-tile files but when you click on it to get the documentation it says:

  This library reads Mapbox Vector Tiles and allows access to the layers
  and features.

As for the Leaflet.TileLayer the documentation says that TileLayer is for rasters.


texts: QGIS 3.10.8 Add features and attributes to exisiting layer using python

texts:
Is there a way to add features from one layer to another while matching attribute names using Python?
Using QGIS 3.10.8
I am trying to add features to an existing layer. I have been using the code below to add the features, this code however does not match the attribute names and populates the added features in the order of the layer I am obtaining features from.
For example: If I have N-1 attributes in the layer I am obtaining the features from, the data will add features, but will add "NULL" to the last attribute columns, while the attributes with "NULL" should have been somewhere else, illustrated in the picture at the bottom.
Is there a way to ensure that features get added to the correct attribute?
layer = QgsVectorLayer(Out, "Out","ogr")
feats = []
for feat in Inn.getFeatures():
    feats.append(feat)
layer.startEditing()
dp = layer.dataProvider()
dp.addFeatures(feats)
layer.commitChanges()


texts:
I have the code below
import arcpy, os

from arcpy.sa import *

arcpy.env.overwriteOutput = True

fc = r'E:\my saves\finaldata.dbf\finaldataD'
# Point feature class, needs to be in a file geodatabase (for the where clause to work)
output_folder = r'C:\GIS\data\testdata'  # Where the rasters will be saved

date_field = 'date'  # field holding the dates (date type, not string)
value_fields = ['CO', 'O3', 'NO2']  # field with the values to interpolate

all_dates = list({i[0] for i in arcpy.da.SearchCursor(fc, date_field)})  # List all unique dates

for date in all_dates:
    # Create a feature layer of each unique date
    where = "{0}=date '{1}'".format(arcpy.AddFieldDelimiters(datasource=fc, field=date_field),
                                    date.strftime('%Y-%m-%d'))
    arcpy.MakeFeatureLayer_management(in_features=fc, out_layer='lyr', where_clause=where)

    # For each field interpolate. This is untested since I dont have spatial analyst
    for fieldname in value_fields:
        outIDW = Idw(in_point_features='lyr', z_field=fieldname,
                     cell_size=10)  # , {power}, {search_radius}, {in_barrier_polyline_features})
        outIDW.save(os.path.join(output_folder, '{0}_{1}.tif'.format(fieldname, date.strftime('%Y%m%d'))))

but in the result I got runtime error
Traceback (most recent call last):
  File "C:/Users/USER/AppData/Roaming/JetBrains/PyCharm2020.1/scratches/arcgis.py", line 15, in <module>
    all_dates = list({i[0] for i in arcpy.da.SearchCursor(fc, date_field)})  # List all unique dates
RuntimeError: cannot open 'E:\my saves\finaldata.dbf\finaldataD'

can anyone tell me how to solve this problem?


texts: How to automate geoserver installation using Chef?

texts:
I need to automate a GeoServer installation on CENTOS 6.x using Chef.

I've found only this question about this issue,  Automating removal of all demo layers from Geoserver, but it's quite outdated

Any suggestions / examples updated?


texts:
Before, to load my .geojson, I had a php file who generated the code :

var vector_limite_comm_white = new ol.layer.Vector({
    source: new ol.source.GeoJSON({
        url: 'geojson1.php',
        defaultProjection :'EPSG:4326', projection: 'EPSG:3857'
    }),
    name: 'NAME 1',
    style: style_white()
})

After, I tried to directly load my geojson file, so I wrote :

var vector_name1    = new ol.layer.Vector({
    source: new ol.source.Vector ({
        url: 'data/json/file1.geojson',
        defaultProjection :'EPSG:4323', projection: 'EPSG:3857'

    }),
    name: 'NAME 1',
    style: style_white()
})

In firebug, at onload event, I saw my geojson file was loaded correctly but, it doesn't appear in my map. So I tried an other method like that :

var test = 'data/json/file1.geojson';
var vectorSource = new ol.source.Vector({
    features: (new ol.format.GeoJSON()).readFeature(test)
});
var vector_overlay = new ol.layer.Vector({
    source: vectorSource,
    style: new ol.style.Style(style_white()
});

It's not a good method too.

Do you have some explications?


texts: How can have access to the inner features (neighbors, edges) of a triagulation remove points from the mesh under PostsGIS

texts:
I've a calculated a triangulation of a point set grouped in certain distance (buffer). For areal classification purposes, I need to calculate:

(a) the points in the unified buffer, which is easy

(b) and need to remove points from the mesh if they are isolated, which is difficult for
    me with postgis board utilities.

I try to find a function which calculates the edge length  for neighbored points (delauney triangulation) and remove points/triangles from the triangulation with large edges.

Details:

Assuming I've a group of points in a plane. The points should be grouped by a set of rules.

1. Collection step Each point with an  overlapping buffer of 75m is belonging to a group.

 -- 1. Add a geometry buffer for each point, here 75 m
 SELECT AddGeometryColumn ('public','stones','b75',32632,'POLYGON',2);
 UPDATE stones set b75=ST_SETSRID(ST_BUFFER(geom,75),32632);

 -- 2. Add a column rid for a region identifier
 ALTER  TABLE stones ADD COLUMN rid INT;

 -- 3. Create a region table for the overlapping
 -- objects of the points
 SELECT * INTO stones_bu75 FROM
  (SELECT
     (A.D).path[1] AS id,    -- region identifier
     (A.D).geom    AS geom   -- unified point buffers
   FROM
     (SELECT st_dump(st_union(bu75)) AS D FROM stones) AS A ) AS S;

  -- 4. update the table of points (stones) and mark/collect each
  --    stone into the region group
  UPDATE stones
    SET rid = C.RID FROM  -- set the region identifier
    ( SELECT
         S.id AS SID ,
         R.id AS RID
      FROM
         stones AS S,      -- get the stone position
         stones_bu75 AS R  -- get the region
      WHERE
         -- check if the stone position is in the region
         ST_CONTAINS(st_setsrid(B.geom,32632),
                     st_setsrid(S.geom,32632))
    ) AS C
  WHERE
      id = C.SID;

As a result I've stone positions grouped by a region id

and can count them (or sort sparse regions out).

2. Compacting step: In a second step I want to sort out points that are isolated, which means a group of points is valid if there are more than 20 points group within a edge length of < 50m to the next neighbor.

So I need to calculate the triangulation within the groups and filter out the large edges:

-- 1. Triangulation of the points for each group
SELECT AddGeometryColumn ('public','stones_bu75','t75',32632,'TIN',2);
UPDATE stones_bu75 SET  t75= T.tri FROM
  (SELECT
     rid AS tid,
     ST_SETSRID(ST_DelaunayTriangles(ST_COLLECT(geom),1, 2),32632) as tri
   FROM
    (SELECT rid, array_agg(geom) as geom FROM stones GROUP BY rid) AS G
  ) AS T
WHERE id = tid;

-- 2. Edge length calculation and filtering
geo=# SELECT id, p, ST_length( ST_MakeLine(sp,ep) )
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT id, p,
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
       -- extract the individual linestrings
      (SELECT id, (G.D).path[1] as p, (ST_Dump(ST_Boundary((G.D).geom))).geom
       FROM
        (SELECT id, st_dump(t75) AS D from stones_bu75) G
       ) AS linestrings
    ) AS segments WHERE ST_length( ST_MakeLine(sp,ep)) < 50;

Question:

Unfortunately, there is no  back references to the id of the points spanning the triangulation (only coordinates). Are there possibilities:

(a) to have a logical access (id) to the points and neighbors of each point from the triangulation set,

(b) how do I remove isolated points (points with no connection to the mesh) from a calculated triangulation.

Best regards Huck


texts: Healthy number of tables in PostgreSQL schema?

texts:
I manage a PostgreSQL + PostGIS database with two schemas of 200 and 120 tables.

I wonder if this is correct or I should organize the tables in more and smaller schemas.

So far, we hadn't had any problems of any kind with the database, despite the annoying scroll to get to the desired table in QGIS file explorer.


texts: Converting OSGB easting and northing to geometry on fly in SQL Server 2008 R2?

texts:
I have a table in a database that contains an easting and a northing for broadband relay stations. I am trying to create a view in SQL Server 2008 R2 that I can use in QGIS and I want to create a geometry field on the fly from the easting and northing. An example is Easting 257306 and northing 142708 (stored in separate fields) converted to POINT(257306 142708) as a geometry.

I've used the following code:

geography::STGeomFromText('POINT(' + CONVERT(varchar(20), Easting) + ' ' +
   CONVERT(varchar(20), Northing) + ')', 27700) AS Geometry

but I get an error that states that the specified SRID must match one of the supported ones in sys.spatial_reference_systems and I've looked and 27700 does not appear in that table.

How do I create this on the fly geometry field in my view?