texts: I've been given two .asc files containing a 1m raster data from a LiDAR survey. I only have access to QGIS.2.14.1 (which I'am using for the 1st time only having previous experience with ArcGIS) and I am struggling to import the files due to their size given that they are 2.2 and 6GB respectively. Is there a function to reduce the size/mosaic the file to enable import in QGIS or any other program?


texts: What causes shiny artifacts in Google Maps imagery?

texts:
Do you know what causes artifacts like this one?

Source: http://maps.google.com/maps?f=q&source=s_q&hl=en&geocode=&q=embelgasse+wien&aq=&sll=37.0625,-95.677068&sspn=60.894251,93.251953&ie=UTF8&hq=&hnear=Embelgasse,+Margareten+1050+Wien,+Austria&t=h&ll=48.187123,16.348997&spn=0.001849,0.004667&z=18


texts:
I am currently using the Dijkstra's algorithm to find the shortest path between multiple latitudes and longitudes on QGIS.

Most of my code is from the Dijkstra's documentation here: http://docs.qgis.org/testing/en/docs/pyqgis_developer_cookbook/network_analysis.html#finding-shortest-paths and I have also referred to a similar issue here: QGIS crash when looping through features in network analysis (dijkstra)

At first, I thought it was Dijkstra's algorithm function that was crashing QGIS. This is still happening, however, I have also found another error. When ran, the console shows "AttributeError: 'NoneType' object has no attribute 'findVertex'", implying that the graph returned from graph = builder.graph() is null. Has anyone run into this issue as well, or any suggestions of how I can move forward?

from qgis.utils import iface
from PyQt4.QtCore import *
from PyQt4.QtGui import *
from qgis.core import *
from qgis.gui import *
from qgis.networkanalysis import *

layers = qgis.utils.iface.legendInterface().layers()
for layer in layers:
    if layer.name() == 'tl_2016_44_tract':
        tractLayer = layer

# building a graph
# all roads are treated as two-way
director = QgsLineVectorLayerDirector(tractLayer, -1, '', '', '', 3)
properter = QgsDistanceArcProperter()
director.addProperter(properter)
crs = qgis.utils.iface.mapCanvas().mapRenderer().destinationCrs()
builder = QgsGraphBuilder(crs)

# extracting coordinates
tractLayer.startEditing()
points = {f.id():f for f in tractLayer.getFeatures()}

for key1 in points:
    lat1 = float(points[key1]["INTPTLAT"])
    lon1 = float(points[key1]["INTPTLON"])
    startPoint = QgsPoint(lat1, lon1)
    for key2 in points:
        lat2 = float(points[key2]["INTPTLAT"])
        lon2 = float(points[key2]["INTPTLON"])
        endPoint =  QgsPoint(lat2, lon2)

        tiedPoints = director.makeGraph(builder, [startPoint, endPoint])
        graph = builder.graph()

        tStart = tiedPoints[0]
        tStop = tiedPoints[1]
        idStart = graph.findVertex(tStart)
        idStop = graph.findVertex(tStop)

        (tree, cost) = QgsGraphAnalyzer.dijkstra(graph, idStart, 0)
        print "Travel time is %s hours" % (dcost[endId])


texts: Calculating polygon areas in QGIS generates NULL results

texts:
I'm having a problem using QGIS and can't find a solution ... I searched this site and found a similar question (QGIS $area field calculator always produces NULL -produces-null), but I couldn't solve my problems with the answers there ...
I have a file and every time I try to calculate the area of the polygons, a lot of them returns as NULL. To calculate the areas I do the following:

Right click on layer > Open Attribute Table
Open Field Calculator
I use the parameters below:

When I click OK to calculate, the table generates several regions with a NULL area value.
In the other publication (linked above) they recommended trying to use the "Lines to Polygons" tools to solve possible problems, or "Check Validity"
The "Lines to Polygons" tool I just can't use ... When I enter it, the "Input Layer" field appears empty and I can't select a layer to run it.

The "Check Validity" tool returns all areas as valid, and the 2 layers generated "error output" and "Invalid output" are empty.
Any idea how I can solve this problem?
Sample file (Shapefile I used to get to this errors):
https://we.tl/t-MuEKE7QF94
I'm using QGIS 3.16.3


texts: When calculating temperature from Landsat is it necessary to use TOA?

texts:
I'm attempting to calculate temperature from Landsat 5 and 8 thermal bands. I've already calculated NDVI, where the NDVI bands were corrected using dark object subtraction.

I see that many people use top of atmosphere reflectance to determine temperature. Would it be incorrect for me to use the DOS bands I calculated to determine temperature, or do I need to go back one step and use TOA reflectance?

Specifically, would it be inappropriate to use DOS radiance for the L lambda term instead of TOA radiance? The equation I'm referring to is the last one on the page I linked.

https://landsat.usgs.gov/using-usgs-landsat-8-product


texts: Is there any shapefile (or other format) of Russia's time zones after 2014 changes?

texts:
According to for example Wikipedia, Russia changed its time zone boundaries in 2014, moving areas between existing zones and defining new ones for the Crimea area.

Is there any shapefile or other geometric data available somwhere which is up-to-date enough to reflect these changes? The only shapefile of Russia's time zones that I could find so far was last updated in April 2013 which was before the change.


texts: gdalwarp leaves black boundary of cropped GeoTIFF with cutline

texts:
This is a new question based on a response on the thread  white page border when converting geoPDF to geoTiff with gdal_translate

I am using gdalwarp as described to do a cutline of GeoTIFF images that have been rendered from geospatial PDFs. The neatline in the PDFs is not accurate but I can easily create UTM coordinates and adjust them accordingly to get the GeoTIFF image that I desire.

This is the CSV file. The coordinates are not square

id,WKT
1,"POLYGON ((410290.5552 4025996.01,410434.6684 4039895.26,421656.358 4039786.037,421530.3048 4025886.935,410290.5552 4025996.01))"

This is the command line:

gdalwarp original.tif original.crop.tif -cutline mycutline.csv -crop_to_cutline

The new image is square and the cropped image is canted inside as is expected. The problem is that the outer crop area is showing up as black. When I tile the image and display it on top of a map base layer, it aligns perfectly except that it has a black edge around it. If the edging was transparent then it would be fine, but the when tiled, it definitely retains the black color.

Here is a sample image
http://www.appalachianbiketrails.org/images/resized.tif

Here is the tile example
http://www.appalachianbiketrails.org/tiles/openlayers.html

I have tried adding the -dstalpha flag as mentioned in other posts and shown in the code snippet below, but that gives the error "Cannot compute bounding box of cutline."

gdalwarp original.tif original.crop.tif -cutline mycutline.csv -crop_to_cutline -dstalpha

On a side note, I found that the direction of the coordinates in my CSV file determined whether the crop was inside or outside of the coordinates.


texts: Export Muliple Layer Selection to Single CSV File in ArcGIS WebApp

texts:
I'm trying to set up a web app that has multiple layers that all have the same fields (each layer is a separate year of data but same data type).

I need to be able to select data points from multiple layers and export the selected features to a single csv file where the records are appended to a single set - not joined since the fields all match up.

Right now, the only widgets I've found allow me to select multiple layers in one go but only export each layer separately. Before I set out trying to build my own custom widget since this seems like a pretty "run-of-the-mill" function for a mapping application, I'm wondering if anyone here knows of a standard or custom widget already out there that achieves this.

I'm working on my app through ArcGIS Web App Builder developer edition.


texts: No route from location ... to location ... Error: No solution found by Network Analyst?

texts:
Why does this warning box appear when I click the Solve button in Network Analyst toolbar?

No route from location "Graphic Pick 46" to location "Graphic Pick 61".
Error: No solution found.


texts: Finding all polygons overlapping polygons from another layer using MapInfo?

texts:
I have two polygon layers. One is All_India_Locality layer which has attributes:{Locality_ID,Locality_Name} and the other is All_India_Pincode layer which has attributes:{Pincode,City}

I want to find all the pincodes that are overlapping the localities. So my output should be a table which has attributes:{Locality_ID,Locality_Name,Pincode}.

There's one catch theat pincodes which are just sharing boundaries with the localities should not be selected. (Therefore intersect query won't work)

Example: For the highlighted locality "8781", only Pincode "574105" should be displayed in Pincode column as it overlaps the locality "8781". Pincodes "576107" and "574118" should not be displayed as they are sharing the boundary.
   If localities are overlapped by multiple pincodes, then multiple records of same locality should be present.

I am using mapinfo. One workaround I found for this query is that I split the localities layer with pincodes layer and then updated the Pincode column with pincode values from pincode layer using query "Pincode.obj contains Localities.obj"

Is there any other way?


texts: time settings in NetCDF file get corrupted by my ArcMap-Program

texts:
I downloaded a NetCDF-file from here: http://eca.knmi.nl/download/ensembles/downloadchunks.php (0.25 reg grid, years 1995-2016, tn), but the time is not displayed correctly when viewing it with ArcMap 10.4 on my laptop (advanced licence) or computer (basic licence):
Instead of giving back each day of all the years, it gives me every 8 Seconds of the first three days. It somehow is reading the date wrongly.

The same file worked fine on the PC of a colleague. Does anyone have a guess, what the problem might be?

I might have corrupted the settings when trying to find out how to work with NetCDF-Files, but how can I straighten that out again?


texts: How to use a different projection than EPSG:3857 with CartoDB or Leaflet?

texts:
I have data that is accurately positioned with the NAD 83 MTQ Lambert EPSG:3798 projection.

If I reproject this data in EPSG:3857 (web mercator) the markers are no longer positioned accurately on the tile layer. Is it possible to use CartoDB, or Leaflet, to use a different projection than EPSG:3857?


texts: Oracle spatial - check polygon vertices counter clockwise

texts:
I need to import a huge number of polygons from a third party (as a csv file- array of x,y's).  I have to check whether the polygon vertices are in clockwise or counter clockwise, finally i have to convert it into counter clockwise and store it in the oracle spatial table as sdo geometry.

I am not able to find any oracle spatial java API to check the array of vertices and convert it into counter clockwise.  I tried the same question in oracle spatial forum, but they said no API for this.

do you have ever worked on this oracle spatial java API, any other thoughts.  ( i have found some methods to validate the ring rotation "http://paulbourke.net/geometry/clockwise/index.html", but if it is from oracle it will be good to implement).


texts: Enabling "Place Overlapping labels" in ArcObjects?

texts:
I am currently using ArcObjects and I'm trying to tick the "Place Overlapping labels" option programatically.

I have been using Standard ILabelEngineProperties and the properties of BasicOverposterLayerProperties.

I just cannot find this particular option


texts: Searching Pixel Location (idrisi) on Another GIS

texts:
I'm Searching similar command Pixel Location (like idrisi) on Another GIS, I, need create raster file  with values lat, long by separately...

Update:

I need a raster file where the value of each cell are Latitude and another raster Longitude, (similar centroid of feature point) but on raster.. So this command on IDRISI Pixel location copy rows & col similar raster e.g srtmdem, and create 2 raster with a distribution for each pixel Latitude an Longitude.


texts: Calculating a vector raster intersection in PostGIS

texts:
I am using CartoDB and trying to calculate the total population (in raster format in popcount) within a specified distance to a line (vector data in mytable). The population data is at about 5km resolution. The approach I took was to vectorize the raster data around the line using ST_DumpAsPolygons, then to Sum() up all of these raster values.

SELECT filter,round(Sum((gv).val)) FROM
(
  SELECT provider, ST_DumpAsPolygons(ST_Clip(rast, 1, mytable.the_geom, true)) AS gv
  FROM popcount, mytable
  WHERE ST_Intersects(rast, mytable.the_geom)
) AS foo GROUP BY foo.filter ORDER BY round DESC

The problem with this approach is that I seem to double count populations that are close to several lines.

Is there a better way to calculate a vector-raster intersection, or is there a modification I can make to this query to avoid double counting?


texts: How to delineate these coloured mandals as different shapefiles using QGIS?

texts: The image shown above is everything I have! Is there a way to make separate shapefiles for all the coloured palletts


texts: Setting Method Parameters for ArcGIS Generating Near Table tool?

texts:
I input firms' coordinates in the arcgis and want to calculate the distance of firms to the nearest polyline feature (roads). I am having a problem grasping the unit of near_dist value. All of them are in hundredth or thousandth like .00475. Most likely the value is expressed in the coordinate system. But I would like them to be shown as in the metric system.

ArcGis has a help file that mentions you need to change the method parameter from planar to geodesic and that conversion will give you an near_dist value in meters.

But I am really having trouble to find the method parameter setting. Can anyone walk me through how to find it. I am using ArcMap 10.2 and windows 32-bit.

output:


texts: Default settings for composition and item properties in print composer

texts:
Is it possible to make some default settings or save the settings for composition and item properties?
Im adding a bunch of different maps to print composer. I want the same settings for every map f.eks the same position and size. When Im adding new maps, I have to change all the settings every time.


texts: How to have two OpenLayers maps on the same page

texts:
Specifically, my question is: how to assign the maps to different divs? In my page it kept placing the two maps under each other, although I defined the divs to be horizontally aligned.
So I have been searching for the answer to this question, when I discovered it by mistake, and I think it is not obvious so I thought I might share it.


texts:
When I add a polygon feature it shows on the map. When I add another one the first one disappears. The code for adding is below:

$.postJSON('@Url.Action("GetLayers")', { "planId": planId },
            function(data) {
                var features = [];
                $.each(data, function(i, item) {
                    var pointList = [];
                    $.each(item.Points, function(y, p) {
                        pointList.push(new OpenLayers.Geometry.Point(p.X, p.Y));
                    });
                    var ring = new OpenLayers.Geometry.LinearRing(pointList);
                    var feature = new OpenLayers.Feature.Vector(new OpenLayers.Geometry.Polygon(ring));
                    feature.url = '@Request.Url.AbsoluteUri';
                    feature.fid = item.Id;
                    feature.geometry.id = "Link_" + item.Link;

                    features.push(feature);
                });
                polyLayer.addFeatures(features);
            }
        );

polyLayer is a Vector layer.

Adding point features to another vector layer works perfect.

Is there something that I'm doing wrong?


texts: Putting a Mapbox-GL canvas into an OpenLayers node results in the Mapbox canvas being slim

texts:
I follow this example but my canvas is slim.
I set a breakpoint here.
const mbMap = new mapboxgl.Map({
  // the node has a clientWidth of 674px.
  container: node
})

The canvas is 15px less wide.
<canvas class="mapboxgl-canvas" tabindex="0" aria-label="Map" width="659" height="981" style="width: 659px; height: 981px; position: absolute; display: block;"></canvas>

The result is a white column to the right of the Mapbox-GL canvas.

The only fix seems to resize the node which is the container for the mapbox-gl canvas. This will resize the canvas accordingly.


texts: Change the data type of a field in multiple shapefiles

texts: Is there a way I can change the data type of a field (possibly using the python window). I have about 25 shapefiles which are all named differently, but the name of the fields I want to change from text to double are the same.


texts: How can i generate parcel extracts on demand through a webservice

texts: Ok, i need to generate parcel extract documents on demand. The backend is Arcgis 9.3.1 on SQL server 2008. I would like to build a webservice (C# .Net) on the server that generates these extracts. My experience is writing with arcobjects to extend Arcmap. But now i just want to access the data and generate without this UI.  What are the options i have? Any suggestions would be great. (p.s. it is not necessary to stick to .NET) I only need some advice or direction. I'll try to describe my problem in general. The extracts will be requested over the internet. So the only thing the system should do is returning an jpg or pdf. That's all. The problem i have now is that i don't now how to do it. Any concept or advise would be great. Sorry if this is still to generic. I have to format the extract from a template so a grid and a north arrow is added. I know how to do this with arcobjects in arcmap but i don't know how to use arcobjects from a webservice and i even don't know if this is possible.


texts: Configuration for GeoWebCache in Geoserver 2.4.0

texts:
I have developed a web map application by using Geoserver 2.4.0 and tomcat 7 in Windows XP.

I have tried to configure the GeoWebCache for WMS layers to speed up the service but unfortunately it didn't work for me. I did something wrong in configuration steps. Kindly Help me to achieve this. I have attached some screenshots for reference.

I have created New Gridset for EPSG:404000

I have tested with XXX.XXX.X.XX:8080/geoserver/gwc/demo Link:

My layer setting is:

And My Geoserver Log is:

at org.apache.catalina.valves.AccessLogValve.invoke(AccessLogValve.java:953)
at org.apache.catalina.core.StandardEngineValve.invoke(StandardEngineValve.java:118)
at org.apache.catalina.connector.CoyoteAdapter.service(CoyoteAdapter.java:408)
at org.apache.coyote.http11.AbstractHttp11Processor.process(AbstractHttp11Processor.java:1023)
at org.apache.coyote.AbstractProtocol$AbstractConnectionHandler.process(AbstractProtocol.java:589)
at org.apache.tomcat.util.net.JIoEndpoint$SocketProcessor.run(JIoEndpoint.java:310)
at java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source)
at java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source)
at java.lang.Thread.run(Unknown Source)
2014-01-18 10:30:31,901 ERROR [geowebcache.GeoWebCacheDispatcher] - Problem communicating with GeoServer http://xxx.xxx.x.xx:8080/geoserver/gwc/service/wms
2014-01-18 10:30:31,901 INFO [geoserver.wms] -
Request: getServiceInfo
2014-01-18 10:30:31,901 ERROR [geoserver.ows] -
org.geoserver.platform.ServiceException: No such style: ODC_Raster
at org.geoserver.wms.map.GetMapKvpRequestReader.parseStyles(GetMapKvpRequestReader.java:1386)

My Default Raster Style:

My Custome Raster Style:

Thanks in advance


texts: R: convert TIN to multipolygon with package sf


texts: Most optimal way to create a (spatial) database for big amounts of data?

texts:
I am using QGIS, and on my way to set up a huge spatial database (geographical), where most of the data is points without geometry), along with other types of data (like lines in DXF and areas in SHP).

Until now I have been using single files for all data, but I begin to believe that this will be rather slow (rendering/querying/calculating) when adding more and more data (?). I have read about a few different methods (like this discussion: Database versus Files for Solo Developer). Unfortunately I haven't found any article that clearly states pros and cons when using "single file data import" compared with using a proper database instead.

At the moment, I am the only one going to use this database, but in the future it will probably be more people, and perhaps even we will need to "install" this database on other peoples computers (not sure how, yet).

I'm rather new to using PostGIS, but in my understanding, to use (for example) a PostGIS database I need to save all my data as Shape files, instead of point text files (CSV)? So this is also an extra step that will take some time (and I feel more comfortable and used to making calculations and changing data in a CSV-file using Excel).

My main question is: what the best way to build a database for this purpose?
A sub-question is: what are the pros and cons when comparing "single file data import"  and a proper database (I plan to use a PostGIS database)?

Some clarifications: First of all, I haven't had much education in GIS engineering, so perhaps some of my expressions and technical terms are somewhat wrong...

"Are you just using the csv to view the coordinates as a point?" - Yes, this is what I wrongfully called "single file import" earlier, the point (coordinates) itself is of course geometry as well.

"Really, what's huge?" - Good question, it's in the eye of the beholder I guess :) In this case it will be a database with around 150 000 points separated into many different filter layers. As it is now, each filtered layer reads from the CSV file everytime it loads, and I experience that to be a rather slow process compared to reading from only one database (with shapefiles) instead (correct me if I'm wrong).
Other than that there is only about 20 000 - 30 000 line/polyline objects (about 200 Mb data).


texts: Very large GeoTIFF dataset to GeoPackage

texts:
I'm trying to mosaic a large GeoTIFF dataset into a GeoPackage with Gdalwarp
(no reprojection). There are 512 .tif files, each with a single 8-bit grayscale
band, and 1 meter resolution. I'm using -cutline to clip the dataset to a
slightly smaller extent.

I'm using the command

gdalwarp --config GDAL_CACHEMAX 35% -multi -of GPKG -wo SKIP_NOSOURCE=YES -co QUALITY=70 -cutline worldbounds.shp -csql "SELECT * FROM worldbounds WHERE country='TheCountry'" *.tif output.gpkg

GdalWarp initially outputs the dimensions as 1 824 323P x 1 662 148L

I expect this operation to take quite some time, but it's now been 7 days and there has been no other output than the message about the dimensions.
GdalWarp.exe has been using 2-3% CPU and 5 GB of memory the entire time. The output.gpkg is still only 80 KB and has not been modified since I started the process. Since GdalWarp usually outputs each filename as it's being processed, I'm
suspecting that something has gone wrong.

From experience with similar jobs, it seems that when the output image is larger than a certain size, the process will not move forward. Does anyone have experience with creating very large mosaic GeoPackages in GDAL (performance hints etc.)? I'd be very interested to hear if you've managed this with other software as well.


texts:
perfect result Shapefiles

How I can convert this closed line to polygon with PostGIS?

I tried many functions like ST_Polygonize, ST_MakePolygon, ST_BuildArea. I got some results, but not a perfect one.

The ArcGIS function Feature to Polygon generated a perfect result.

The shapefile is available on request.


texts: Can't export customer's shapefile to latitude and longitude

texts:
I'm very new to this application so bear with me.  I hope that I am missing something very elementary.

I have a shapefile supplied by my customer.  I need to get the latitude and longitude for each of the points.  From the export it appears there are 140 points.

I looked at the attribute table and only the OBJECTID, Shape_Leng are defined.

I have tried to do the steps in this post but I only get a file that looks like:

OBJECTID,Shape_Leng
1,17715.43075900000
2,20656.19732000000

What I was hoping to get a a file that looks like:

OBJECTID, latitude, longitude

The customer sent me a zip file with all of this data:

customer_x.dbf
customer_x.qix
customer_x.shp.xml
customer_x.prj
customer_x.shp
customer_x.shx

I have been opening the .shp file and haven't been able to open any of the other files with QGIS.


texts: Standardized Precipitation Index with gamma function

texts:
Inspired from Climate Engine platform, I want to calculate Standardize Precipitation Index based on this paper, where each rainfall datasets are fitted to the gamma function to define the relationship of probability to rainfall.

I have try to write a code in GEE to calculate 3-month SPI, so far using simple equation where SPI = (current rainfall - LongTermAverage) / StDev spi = rainfall_now.subtract(ltamean).divide(ltastd). I have difficulties on how to fitted the data to gamma function in GEE, is there any reference or example about this?

Map.setCenter(109.738, -7.347, 7);

// We will compare for a 1-month period of CHIRPS data
var start = 31+28;

// Long term average (LTA) can be generated directly from the pre-calculated Rainfall
var lta = ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD')
  .filter(ee.Filter.dayOfYear(start, start + 31))
  .filterDate('1981-01-01', '2019-12-31');

var ltamean = lta.select('precipitation').mean();
var ltastd = lta.reduce(ee.Reducer.stdDev());

// same for this year
var latest = ee.ImageCollection('UCSB-CHG/CHIRPS/PENTAD')
  .filter(ee.Filter.dayOfYear(start, start + 31))
  .filterDate('2020-01-01', '2020-05-01');

var rainfall_now = latest.select('precipitation').mean();

// add the various layers and layer combination
Map.addLayer(rainfall_now, rainfall_viz, 'JFM 2020', false);
Map.addLayer(ltamean, rainfall_viz, 'LTA', false);
Map.addLayer(ltastd, rainfall_viz, 'STDev', false);

// Calculating rainfall anomaly and standardized precipitation index
var anomaly = rainfall_now.divide(ltamean).multiply(100);
var spi = rainfall_now.subtract(ltamean).divide(ltastd);

// Define a colour table
var red2blue = ['a60027','de3e2c','fa8b50','ffd485','e1e1e1','d3edf5','90c4de','5183bd','313694']
var rainfall_viz = {min:0, max:25, palette: ['B8810B', 'CB8C08', 'FFE47E','FFF2C1', 'DAF7A6', 'D2FFFF', '2DE7E7', '00B9D6', '008BD6', '006399']};
var rainanom_viz = {min:0, max:200, palette: red2blue};
var spicolor = ['760005','ec0013','ffa938','fdd28a','fefe53','ffffff','a2fd6e','00b44a','008180','2a23eb','a21fec'];
var spi_viz = {min:-2, max:2, palette: spicolor};

// Add the image to the map using both the color ramp and interval schemes.
Map.addLayer(anomaly, rainanom_viz, 'Anomaly (JFM2020-LTA)');
Map.addLayer(spi, spi_viz, 'SPI JFM2020');

For reference, attached is the spreadsheet file, that I frequently used to calculate SPI3 that incorporate gamma function and normal inverse in the calculation.


texts: OpenLayers measurment path in 31700 projection

texts:
I have this code who works great for distance measurment on 900913 projection with 4326  display projection (bing maps + layer WMS):

measureDistanceControl = new OpenLayers.Control.Measure(
                  OpenLayers.Handler.Path, {
                    persist: true, // 1
                    immediate: true, // 2
                    displaySystem: 'metric', // 3
                    eventListeners: {
                        'measure': measure2,
                      'measurepartial': updateDisplay // 4
                    }
                  }
                );

                function updateDisplay( event ) { // 5
                    var mesaj= event.measure + ' ' + event.units;
                    document.getElementById('output').innerHTML=mesaj;
                };

                function measure2(event)
                {
                var mesaj= event.measure + ' ' + event.units;
                    document.getElementById('output').innerHTML=mesaj;
                }

    }

    function measure()
    {
        map.addControl( measureDistanceControl );
        measureDistanceControl.activate();
        document.getElementById("cancel").style.display = "inline";
    }
    function cancel_measure()
    {
        measureDistanceControl.deactivate();
        map.removeControl( measureDistanceControl );
        document.getElementById("cancel").style.display = "none";
    }

When I insert the code to a page with a layer and 31700 projection it sends NaN error. Code:

 map = new OpenLayers.Map({
    div: "map",allOverlays: true, projection: new OpenLayers.Projection("EPSG:31700"),
   allOverlays: true,maxResolution: 0.4651995555350368, maxExtent: bounds,
    layers: [

        data=new OpenLayers.Layer.WMS("data","http://127.0.0.1:8080/geoserver/data/wms",{layers: 'morino:<?php echo $data; ?>_space',transparent: true, format: 'image/png', <?php if ($tic!="") { ?>CQL_FILTER: "id_ddd=<?php echo $eras; ?>"<?php } ?>},{isBaseLayer: true, opacity: 0.9, singleTile: true }),
        wfs

    ],
    controls: [
        new OpenLayers.Control.Navigation({
            dragPanOptions: {
                enableKinetic: true
            }
        }),
        new OpenLayers.Control.PanZoom(),
        new OpenLayers.Control.Attribution()
    ],
    //center: new OpenLayers.LonLat(26.1473586305531, 44.4431676053183).transform(new OpenLayers.Projection("EPSG:4326"), new OpenLayers.Projection("EPSG:31700")),
    zoom: 15
});

I've searched on Google but I didn't find anything to help me. I also read OpenLayers Cookbook - but the examples are only for 4326 projection. Or maybe OpenLayers doesn't support 31700 projection (or other then 4326)?

EDIT Solution found. Add this in header

<script type="text/javascript" src="js/proj4js/proj4js-compressed.js"></script>
    <script type="text/javascript">
     Proj4js.defs["EPSG:31700"] = "+proj=sterea +lat_0=46 +lon_0=25 +k=0.99975 +x_0=500000 +y_0=500000 +ellps=krass +towgs84=33.4,-146.6,-76.3,-0.359,-0.053,0.844,-0.84 +units=m +no_defs";
    </script>

Use proj4js v1 because proj4js v2 has a problem with openlayers2. I used http://svn.osgeo.org/metacrs/proj4js/trunk/lib/proj4js-compressed.js


texts: Google Maps/Google Earth - KML Overlays and Custom Orientation

texts:
I am am trying to get the floor plans for our various buildings imported into 4 separate KML files, one for each floor level so that I can toggle them with javascript. I'm able to overlay our floor plan images and properly scale them in Google earth, but once exported the KML doesn't display properly.

Please take note that whatever the orientation of the Google map is, I can't get my overlays to display properly. On the actual API based map, I have it set to be oriented with West at the top instead of North, which seemingly throws us another curve ball, as when I attempt to load the KML per the API's docs and refresh the API based map, I end up with nothing but "There are no images for this zoom level."

My KML:

<?xml version="1.0" encoding="UTF-8"?>
<kml xmlns="http://www.opengis.net/kml/2.2" xmlns:gx="http://www.google.com/kml/ext/2.2"     xmlns:kml="http://www.opengis.net/kml/2.2" xmlns:atom="http://www.w3.org/2005/Atom">
<GroundOverlay>
<name>Untitled Image Overlay</name>
<Icon>
    <href>http://img694.imageshack.us/img694/2600/walther.png</href>
    <viewBoundScale>0.75</viewBoundScale>
</Icon>
<LatLonBox>
    <north>37.21924494476235</north>
    <south>37.21876115002708</south>
    <east>-93.2637540015874</east>
    <west>-93.26445522810076</west>
    <rotation>87.56327787656277</rotation>
</LatLonBox>

My API v3 Javascript:

<script type="text/javascript">
    function initialize() {
        var options = {
            center: new google.maps.LatLng(37.22090,-93.26371),
            zoom: 18,
            minZoom: 18,
            maxZoom: 21,
            heading:-90,
            mapTypeId: google.maps.MapTypeId.SATELLITE
        };
        var map = new google.maps.Map(document.getElementById("map_canvas"), options);

        //here's where it breaks!!!
        //it just won't load the kml...or it is and just breaks beyond reason!
        var ctaLayer = new google.maps.KmlLayer('test.kml');
        ctaLayer.setMap(map);
    }
</script>

Here's what it looks like in Google earth (correct):

http://img811.imageshack.us/img811/6348/gearth.jpg

Here's what it looks like in Google maps (for testing only, using it to load the KML...until I also figure out how to get my API based map to even display the layer without breaking):

http://img838.imageshack.us/img838/2237/gmapsr.jpg

Does anyone here have any idea what I need to do to get the display working properly, with regular google maps and more importantly my API based google map?


texts: Add legend from raster catalog? ArcGis 10.1

texts:
Does anyone have a solution on how to add the legend from a raster catalog?

It does not show up, either in version 10.1 or 10.2.1. It only shows the name "Legend".


texts: How to replace a space with an underscore when initiating image export in Google Earth Engine

texts:
I am running several thousand polygons through a Google Earth Engine code and many of them have names with spaces in them and when I try to run the initiate image export I get "Task name may contain alphanumeric characters, underscores, and dashes. No spaces." It would be impractical for me to go through and add and underscore to them in the attribute table.

So in this one case I am trying to run a fire with the name Flying R. In the attribute table it has a space in the name so it will not successfully run. Is there a way I can code it so that it will automatically insert an "_" wherever there is a " "?

This is the code I though would be relevant. Let me know if you need more.

var fires = ee.FeatureCollection("users/mtd25/fires/AZ3143311069520170614_FLYING_R");
.
.
var fireID    = ee.List(fires.aggregate_array('Fire_ID')).getInfo();
print(fireID)
var nFires = fireID.length;
var Fire_Name = ee.List(fires.aggregate_array('Fire_Name')).getInfo(); //For Megan's dataset
.
.
for (var j = 0; j < nFires; j++){
  var id   = fireID[j];
  var Name = id;
  var fireExport = ee.Image(indices.filterMetadata('fireID', 'equals', id).first());
  var fireBounds = ee.Feature(fires.filterMetadata('Fire_ID', 'equals', id).first()).geometry().bounds();

  for (var i = 0; i < nBands; i++) {
    var bandExport = bandList[i];
    var exportImg = fireExport.select(bandExport);
    Export.image.toDrive({
      image: exportImg,
      description: Name + '_' + Fire_Name + '_' + bandExport,
      fileNamePrefix: Name + '_' + Fire_Name + '_' + bandExport,
      maxPixels: 1e13,
      scale: 30,
      crs: "EPSG:4326",
      folder: 'fires',
      region: fireBounds


texts: Getting Landsat 7 images scanline errors in ArcGIS Desktop?

texts:
I have used "Fix Landsat 7 Scanline Errors" tool to remove the stripes from the raster. Though the black scan lines are removed but the raster gets scratches. See
 screen shot below.

How can I solve this or is there any other way to remove the scan lines?


texts: How to put image over the layer in OpenLayers

texts: How to set the Latitude and Longitude range that focuses on a specific region in Mercator Projection when the layer loads in OpenLayers API ? Is there any way to put image over the layer?


texts: Mapping roadway distance to polygon using ArcMap?

texts:
I need to map the roadway distance from a point to a polygon (or polygon to polygon).

Is this possible in ArcMap?

I know I can't add polygons as a location in Network Analyst so I'm not sure if it's possible.


texts: Copying and pasting features using PyQGIS?

texts:
I have created a table join between two shapefiles which I've named 'OSMM' and 'Mastermap'. I want to copy and paste all features (polygons) without a match from the target layer (Mastermap) into the parent shapefile (OSMM). I am relativley new to PyQGIS and can't find any up-to-update code for QGIS 3.6.

Here is my code so far.

#1. loads the SHP files
uri = 'C:/Users/xx/Downloads/OSMM.shp'
vlayer = iface.addVectorLayer(uri,"","ogr")
uri = 'C:/Users/xx/Downloads/Mastermap.shp'
vlayer = iface.addVectorLayer(uri,"","ogr")

#2. Joins "InvoiceNum" from OSMM to Mastermap
layerToJoin = QgsProject().instance().mapLayersByName('OSMM')[0]
target = QgsProject().instance().mapLayersByName('Mastermap')[0]
myJoin = QgsVectorLayerJoinInfo()
myJoin.setJoinFieldName('InvoiceNum')
myJoin.setTargetFieldName('InvoiceNum')
myJoin.setJoinLayerId(layerToJoin.id())
myJoin.setUsingMemoryCache(True)
myJoin.setJoinLayer(layerToJoin)
myJoin.setJoinFieldNamesSubset(['InvoiceNum'])
target.addJoin(myJoin)

#3. Selects NULL features from Mastermap
params = { 'EXPRESSION' : '\"OSMM_InvoiceNum\" is NULL', 'INPUT' : 'C:/Users/xx/Downloads/Mastermap.shp', 'METHOD' : 0 }
processing.run("qgis:selectbyexpression", params)

So how do I perform the PyQGIS version of Ctrl+C on 'Mastermap' and Ctrl+V into 'OSMM'?

Presumably I will then need to use the following code:

#5. Save the changes
OSMM.commitChanges()


texts: i am carrying out a research on dust, i have a low programming skills, but in order for me to do this effectively, i need to determine the soil moisture of the study area. But the soil moisture files i have are in GeoTiff format, can anyone suggest how to convert these files from a GeoTiff to NetCDF file.


texts: Splitting a line into multi-line at desired points in arcgis 10.1

texts:
I am currently working on a project in ArcGIS 10.1 Advanced license regarding the a single facility multiple destination analysis. But, I encountered a problem. I have to separate line features into multiple segments at desired points. And also I wanted each segment identified individually after the separation of lines at points. I have to do this for network analysis.

The lines and points are shapefiles. The points are generated at using fixed distance tool. To evaluate time for each segment between the points, I need distance between those points. So, please suggest me any method to achieve this.


texts: Counting points in polygon containing holes using QGIS?

texts: Basically, I have these area (part of a bigger map) containing many holes (e.g. lakes) I have thrown on the map some 1000 random points and I want now to count those that fall in that area. I already tried with "Data Management tools>Join by attributes" but I receive a geometry error. I also tried the "Vector>Count Points in a Polygon" but I still receive the geometry error. As I said before, I guess the issue is because the area I want to overlap with the random points contains holes. I would be really surprised to know that that one is the only reason why QGIS cannot count the points but I have to say I am a new to that software.


texts:
I have a road network where I have to cover every road for visual inspection. I need to design a route map for this task.

What tools are available to a QGIS user? How would you describe the workflow in this situation?


texts:
I've had a request recently to generate a map which includes UK Postcode districts, ie aggregations of postcodes based on the first section of the code (NW2, B1, E11 etc.).

From what I've read, these aren't available freely (Postcode Shape Files for the UK).

Is there a way of scripting aggregation of boundaries in either QGIS (or possibly R)? I've done this sort of aggregation before using MapBasic but no longer use MapInfo (plus it gave really messy boundaries).

I'm thinking I may use Voronoi regions rather than the actual postcode boundaries as this will be easier in terms of file handling (my postcode shapefiles are split across lots of files so easier just to run Voronoi on the OS Codepoint data set, plus I imagine this will give cleaner, if less accurate, boundaries)


texts: Field mappings object not functioning with Merge_management

texts:
I'm trying to limit the number of fields when I perform a merge to two. Below is the script I have created to do so. However, the merged feature class carries over all fields, despite my use of a fieldmappings object.

di is a dictionary with a shape type ("polyline", "polygon", etc) as its key, and a python list of feature classes as its value. Any ideas?

relevant code:

fms = FieldMappings()
jnFm = FieldMap()
sourceFm = FieldMap()
jnFm.addInputField(fc, jobField)
sourceFm.addInputField(outFcFullPath, "FC_SOURCE")
fms.addFieldMap (jnFm)
fms.addFieldMap (sourceFm)

print "\nMerging"
for shapeType in di:
    featureClasses = di[shapeType]
    mergeFC = os.path.join (outGDB, "All_" + shapeType)
    Merge_management (featureClasses, mergeFC, fms)
    print "Created: All_" + shapeType


texts: Unable to load DLL 'ArcGISVersion.dll': The specified module could not be found. (Exception from HRESULT: 0x8007007E)

texts:
I recently migrated my desktop from ArcGIS 10.1 to 10.3.1 and visual studio from 2010 to 2013. I am now getting a problem trying to run an Engine application. It builds fine but when I try to debug it I get the above error. The error occurs right at the start when trying to get a list of installed environments "For Each item As ESRI.ArcGIS.RuntimeInfo In irs":

   Private Sub MyApplication_Startup(ByVal sender As Object, ByVal e As Microsoft.VisualBasic.ApplicationServices.StartupEventArgs) Handles Me.Startup

        Dim irs As Generic.IEnumerable(Of ESRI.ArcGIS.RuntimeInfo)

        irs = ESRI.ArcGIS.RuntimeManager.InstalledRuntimes
        For Each item As ESRI.ArcGIS.RuntimeInfo In irs
            System.Diagnostics.Debug.Print(item.Path)
            System.Diagnostics.Debug.Print(item.Version)
            System.Diagnostics.Debug.Print(item.Product.ToString())
        Next
        Dim bInit As Boolean
        bInit = ESRI.ArcGIS.RuntimeManager.Bind(ESRI.ArcGIS.ProductCode.EngineOrDesktop)

    End Sub

From what I've read the problem is not that the ArcGISVersion.dll cannot be found (it imports into the solution no problem) it is that this .dll seems to call another .dll that cannot be found. The stack trace from the error is:

       at ESRI.ArcGIS.NativeInterop.DllGetClassObject(Guid& rclsid, Guid& riid, Object& ppv)

   at ESRI.ArcGIS.RuntimeManager.GetArcGISVersion()

   at ESRI.ArcGIS.RuntimeManager.<get_InstalledRuntimes>d__0.MoveNext()

   at SanGIS_RoadMaintenance.My.MyApplication.MyApplication_Startup(Objectsender, StartupEventArgs e)
in R:\Work\Sangis\SanGISEditTools\SanGIS_RoadMaintenance\SanGIS_RoadMaintenance\Utilities\ApplicationEvents.vb:line 24

I'm not familiar with ESRI.ArcGIS.NativeInterop and I don't know how to add a reference to it.

Also I created a new ArcEngine project from scratch and that runs fine, no errors when debugging or binding to the license manager.

Does anybody know why this app might be failing after migration?


texts:
I have table network which has about 2 million records (only 3 countries in Europe). When I run Shorthest-path it is very slow, it takes 3 min to take records only.
I tried with dijkstra-sp-delta-directed but it is also slow, about 2-2.5 min.
This sp is using bounding box and gist index so it should be faster but it is not.
In the future I should add some other country so table will be very large.

What could I do to make this routing faster? As far as I know Postgre doesn't have option to save temp table in RAM, it works only in one session.

I have 4 GB RAM, I work with OSM and Postgre because it is free, I realy don't know is there any solution for this problem? I am not in position to buy some expensive solution, I have demand to work with open sources. Of course everyone expects quick and cheap solution that works like Google :)

Update:

For the fastest route (avoid toll , avoid highway) my query is:

SELECT * FROM shortest_path(
   'SELECT gid AS id,
           source::int4 AS source,
           target::int4 AS target,
           time + cost::float8 AS cost,
           time + reverse_cost::float8 AS reverse_cost
    from network where toll is null and road_class!=1 ',
    257027, 276521, true, true)

For the shothest route (avoid toll , avoid highway) my query is:

SELECT * FROM shortest_path(
   'SELECT gid AS id,
           source::int4 AS source,
           target::int4 AS target,
           length + cost::float8 AS cost,
           length + reverse_cost::float8 AS reverse_cost
    from network where toll is null and road_class!=1 ',
    257027, 276521, true, true)


texts: How to use ArcGIS FeatureToLine with Managed Assembly

texts:
I am writing an ArcMap Add-In to combine multiple shapefiles and then run some custom python scripts on the resulting shapefile. However I cannot seem to get the Geoprocessor working to run the FeatureToLine tool.  I am using the ArcGIS SDK 10.4 for .NET.

private void Button5_Click(object sender, EventArgs e)
    {
        Geoprocessor gp = new Geoprocessor();
        gp.OverwriteOutput = true;

        try
        {
            FeatureToLine FTLTool = new FeatureToLine();

            List<string> inputs = new List<string>();

            string[] ins = new string[] { @"D:\\Code projects\\Network-WZRD\\Network Wzrd\\Data\\Centerlines\\Roadway\\Delaware.shp",
                @"D:\\Code projects\\Network-WZRD\\Network Wzrd\\Data\\Centerlines\\Trails\\Trails.shp" };

            FTLTool.in_features = ins;
            FTLTool.out_feature_class = @"D:\\Code projects\\Network-WZRD\\Network Wzrd\\atestout.shp";
            ITrackCancel pTrackCancel = new CancelTracker();
            pTrackCancel.CancelOnKeyPress = true;

            gp.Execute(FTLTool, pTrackCancel);
        }
        catch(Exception ex)
        {
            MessageBox.Show("Error\nAn error occured while attempting to merge shapefiles.\n" + ex.Message);
        }

    }//button5_Click -- Final Load

I know I am incorrectly inputing my shapefiles.  But I couldn't find any examples of tools that used multiple inputs in the In_Feature field. I get a general error when it gets to execute the tool

Error HRESULT E_FAIL has been returned from a call to a COM component

The questions I have are: How can I input multiple shapefiles into the tool? Does the output need to be an existing shapefile object within the code, or can it just be a path to where I would like the resulting shapefile to be put?

EDIT:
I tried Hornbydd's suggestion, and used the exact example in the sample code and got the same error code.

    public void ExampleGPUsingMultipleInputs()
    {

        // Initialize the geoprocessor.
        Geoprocessor GP = new Geoprocessor();

        // Initialize the Union tool.
        Union uniontool = new Union();

        // List all feature classes in the workspace.
        GP.SetEnvironmentValue("workspace", @"D:\Code projects\Network-WZRD\Network Wzrd\Data\Centerlines\Trail");
        IGpEnumList fcColl = GP.ListFeatureClasses("*", "", "");

        // Create a StringBuilder object and add the feature classes as a semicolon-delimited list.
        StringBuilder sb = new StringBuilder();
        string inputfeatures = fcColl.Next();

        while (inputfeatures != "")
        {
            System.Diagnostics.Debug.WriteLine(inputfeatures);
            sb.AppendFormat("{0};", inputfeatures);
            inputfeatures = fcColl.Next();
        }

        // Execute the Union tool with multiple inputs creating landuse.shp.
        uniontool.in_features = sb.ToString();
        uniontool.out_feature_class = "landuse.shp";
        GP.Execute(uniontool, null);

    }


texts: Is there any polygon select map tool for QgsMapCanvas?

texts: I am trying to create a plugin that has a QgsMapCanvas of itself. The question is that I can't seem to find the appropriate library to use polygon select feature in the QgsMapCanvas. I know that iface has actionSelect().trigger() function to implement the required but that allows selection only over the QGIS window and not on the plugin's map canvas.


texts: Manually add a Z variable in Qgis onto a vector layer?

texts:
I am trying to manually add one constant Z-variable (e.g. 200meters) across an entire Vector Layer in Qgis so that I may eventually use the Interpolation plugin to generate a Tin.

I have not been able to find any information on how to do this on the website. Is this possible?


texts: How do you test for a NULL string  field in the Field Calculator with Python?

texts:
I am trying to perform an operation on a string field with a condition that it contains data. Referencing None with Python does not appear to work when the field is nullable and contains <Null> for its data. Am I missing something? I have tried len() and "<Null>" and "Null" as well. What do?

The code, for kicks:

def apt(addr,aptnum):
    if aptnum != None:
        a = addr.split(aptnum,1)[-1]
        return a
    else:
        return addr

Update: It would appear that my problem is that field calculations are not processing on any rows where a null column is tested. That is outside the scope of this question, however, so I will leave it for the time being.


texts: How to calculate protrusion in raster data using ArcGIS?

texts:
I have two raster data set.

Raster 1 is soil depth data (10*10 size )
Raster 2 is ground surface data.(10*10 Cell size)

In most cases, Soildepth data is beneath the ground but in some cases, soil depth data is protruded and it's above the ground surface. Using Arcgis I want in the result Raster 3, soil depth data value and value of ground surface data wherever protrusion is present?

The protrusion part is highlighted in Red closed to the area.


texts:
I am a novice programmer, but I know enough to be dangerous. Here's my issue and hopefully someone knows the answer.

I am developing an application with 4 tabs on the left and a map from ArcGIS Online(AGOL) on the right. What I want to do is when tab1 is active (focused), I want the webmap to be webmap1. When the user clicks tab2 I want the webmap to change to webmap2, and so on.

This is seemingly a simple task, but I cannot find enough documentation to understand how to make this happen. Here's a bit of code starting the initial map that works correctly:

var webmap =  "e1cff258e4674ab1985e2303106b0cc8";
arcgisUtils.createMap(webmap,"mapDiv").then(function (response){
map = response.map;
});

What if want is to make myvariable webmap dependent upon the active tab:

var webmap = (if tab1.active){"e1cff258e4674ab1985e2303106b0cc8"}
             (else if tab2.active){"3dea9137545e41eb8d4877ac8970294f"} etc. etc. etc....

Does someone have any idea how to do this?


texts: ArcMap: interpolation for single polygons individually

texts: The red polygons in the image represent the boundary of infills. The raster (colored area) within the polygons is interpolated from the point shape file (black dots). The image shows an example where interpolation is performed across polygon boundaries. However, the values within polygon A should not affect the interpolation in polygons B-E and vice versa. Is there a straightforward way to restrict the interpolation to the polygons? I'd like to avoid to disassemble the points for each polygon and to assemble the result afterwards.


texts: How to populate gdalwarp with a dynamic value in QGIS geoprocessing?

texts:
I am attempting to call gdalwarp (python binding gdal.Warp()) as a script inside QGIS modeler (QGIS 2.18.10, gdal 2.2.0).  The goal is to correlate any loaded raster with a UTM zone grid, get the correct EPSG code from that grid, and pipe that code as a dynamic value into gdal.Warp().

To do this, I have written a geoprocessing script because the existing Warp tool does not allow one to dynamically set the EPSG code.

I have looked through the QGIS geoprocessing documentation, but it is a bit thin on examples, which makes troubleshooting difficult.

Here is what I have so far:

##Raster_Layer=raster # Raster to be reprojected
##Joined_Layer=vector # Layer containing the correct EPSG code under the 'EPSG' field
##Output_Raster=output raster

from qgis.core import *
from osgeo import gdal

def get_epsg(lyr):
    for field in lyr.fields():
        fname = field.name()
        if fname == 'EPSG':
            epsg_code = field.value()
    return epsg_code

layer = processing.getObject(Joined_Layer)
srs = get_epsg(layer)

rlayer = processing.getObject(Raster_Layer)
input_raster = gdal.Open(rlayer)
output_raster = Output_Raster
srs_string = 'EPSG:{}'.format(str(srs))
gdal.Warp(output_raster, input_raster, dstSRS=srs_string)

I am sure I'm missing some QGIS-specific wrappers for the api calls, but have no clue where to look for a thorough documentation on using gdal in QGIS geoprocessing scripts (if such a thing exists).


texts: CSS trick to add margin on Carto's embedded iframe (for mobile)

texts:
I am using carto on my website through embedded iframes, and I have an issue on mobile : the map takes the full width of the screen, making it hard to scroll to the bottom of the page, as can be seen on the pic.

Could anyone give me a CSS workaround to give all the iframes a left and right margin? Given that I cannot add the margin to the containers directly.

<iframe width="100%" height="520" frameborder="0"    src="https://topraf.cartodb.com/viz/5c19c6fe-34bc-11e6-a91c-0e31c9be1b51/embed_map" allowfullscreen webkitallowfullscreen mozallowfullscreen oallowfullscreen msallowfullscreen></iframe>


texts: Trying to find my GPS corners on property

texts: Based off the legal description Legal Description: TR IN 12-10-72 DESC: BEG AT NW COR, TH S 89 53' E 857.5 FT, S 0 3' W 1778 FT, N 89 53' W 857.5 FT, N 0 3' E 1778 FT TPOB CONT 35.0008 AC M/L..... I get that the property is basically 1778 feet by 857 feet, but it doesn't give me the four coordinates for the corners ...... How do I find these as I can enter them into my Garmin and come relatively close to the actual corners. I have no reason to know the exact corners as I wouldn't be building anywhere near the property lines I's more for when we are out there and just don't want to stray onto other peoples land.


texts: Seeking boundaries of Canadian Location Code regions used for Meteorological Services?

texts: I have looked into the GADM datasets. The problem is that Canada has a very non-standardized level-2 Administrative Zone identification across the various Territories and Provinces. As a result, their Meteorological Weather Forecasting uses a unique system of boundaries they created. Info on this here. Visually, you can see these boundaries here on a map. I am looking for datasets to reproduce those regions (ideally shapefile format, but I'm open to other formats so long as I can generate/transform the geometries from it using something like QGIS).


texts:
Is there a way to tell if a map service is using an SDE database as its data source, through either the REST or SOAP APIs?

I'm working with an app that queries a map service, and it is not returning the correct results if the service is based off of shape files or personal geodatabases.


texts: Deriving grid of points from gridded vector dataset using R?

texts: I have a large gridded vector dataset (see extract below) and would like to derive a grid of points centered in each cell of the grid. Is there a way to do this in R?


texts: A column/attribute has a different length and precision than in the PDF output of the same table in QGIS

texts:
The fields populated by running my script have a different length and precision in QGIS, than what appears in the output of the attribute table with a PDF. I should also mention that before the output is created the attribute table is joined with another attribute table, in case that makes a difference.
Creating and populating the new attributes.
primary = QgsField('Primary', QVariant.Double, "double",  2,  2)
huff_model.addAttribute(primary)
index_pri = huff_model.fieldNameIndex('Primary')

for feature in huff_model.getFeatures():
    probability = feature.attributes()
    ctuid = feature["CTUID"]
    if not probability[col_num]:
        print ("This mall does not exist")
    else:
        # Primary market values will go into this column
        if probability[col_num] >= .6:
            #print("Primary market area: {} ({})".format(ctuid, probability[col_num]))
            huff_model.changeAttributeValue(feature.id(), index_pri, probability[col_num])
        # Secondary market values will go into this column
        if .4 <= probability[col_num] < .6:
            #print("Secondary market area: {} ({})".format(ctuid, probability[col_num]))
            huff_model.changeAttributeValue(feature.id(), index_sec, probability[col_num])

# Delete all other features (rows) that do not have a probability in the primary and the secondary column
expr = QgsExpression("Primary is NULL and Secondary is NULL")
for f in huff_model.getFeatures(QgsFeatureRequest(expr)):
    huff_model.deleteFeature(f.id())

huff_model.updateExtents()
huff_model.commitChanges()
huff_model.updateFields()

Loading Print Composer and creating PDF
# Create the attribute using the map composer in QGIS and show the output in PDF format.
mapRenderer = iface.mapCanvas().mapRenderer()
c = QgsComposition(mapRenderer)
c.setPlotStyle(QgsComposition.Print)
x, y = 0, 0
w, h = c.paperWidth(), c.paperHeight()
composerMap = QgsComposerMap(c, x, y, w, h)

table = QgsComposerAttributeTable(c)
table.setComposerMap(composerMap)
table.setScale(1)
table.setVectorLayer(huff_model)
table.setMaximumNumberOfFeatures(huff_model.featureCount())
c.addItem(table)

# Checks for the existence of a file with the same name and removes it.
if os.path.isfile(PDF_file):
    os.remove(PDF_file)

# Set the properties of the PDF which will take and show the attribute table
printer = QPrinter()
printer.setOutputFormat(QPrinter.PdfFormat)
printer.setOutputFileName(PDF_file)
printer.setPaperSize(QSizeF(c.paperWidth(), c.paperHeight()), QPrinter.Millimeter)
printer.setFullPage(True)

pdfPainter = QPainter(printer)
paperRectMM = printer.pageRect(QPrinter.Millimeter)
paperRectPixel = printer.pageRect(QPrinter.DevicePixel)
c.render(pdfPainter, paperRectPixel, paperRectMM)
pdfPainter.end()

In the image below you can see that the attribute corresponding to primary looks fine. There appears to be the correct precision and length.
Finished Script Attribute Table

However as seen below, in the PDF output of the attribute table, there are far more decimals than what's seen in the attribute table in QGIS.
Same Attribute Table in the PDF Output

There were no new calculations used to create the numbers in the "Primary" column. The original data had a length of 23 and precision of 15, and I want to reduce it to a length of 2 and a precision of 2. The second image with the PDF confirms this as there are exactly 15 decimals places. My guess is the output is still using the "old" format for some reason, I'm not sure why though. I want to make sure that any changes that I make to any field will be updated in the output, in this particular case, the changes to length and precision.


texts: ERROR gdal 3.0.4 comipation : libtool: error: 'port/cpl_base64.lo' is not a valid libtool object

texts:
I'm trying to compile GDAL 3.0.4 on fresh installed Ububtu 18.04
Previously I compiled and installed Proj
$ proj
Rel. 7.0.1, May 1st, 2020

Then installed FileGDB_API
After that I downloaded gdal 3.0.4 and execute configure
 ./configure --with-python \
--with-geos --with-geotiff --with-jpeg --with-png --with-expat --with-libkml --with-xerces-c \
--with-fgdb=$HOME/src/FileGDB_API-64gcc51 \
--with-openjpeg --with-pg \
--with-curl --with-spatialite

make

Them I executed
cd gdal-3.0.4; make install
Error here:
libtool:   error: 'port/cpl_base64.lo' is not a valid libtool object
GNUmakefile:66: recipe for target 'libgdal.la' failed
make[1]: *** [libgdal.la] Error 1
make[1]: Leaving directory '/home/user/src/gdal-3.0.4'
GNUmakefile:77: recipe for target 'check-lib' failed
make: *** [check-lib] Error 2


texts: Remove MapServer comments in GetCapabilities request

texts:
MapServer adds a XML-comment right before the <Service> tag in a GetCapabilities repsonse like this:

<!-- MapServer version 7.4.2 OUTPUT=PNG OUTPUT=JPEG [...] -->

Is it possible to remove this comment somehow?


texts: ArcSDESQLExecute is not able to DELETE rows in one database but works in other

texts:
I have two ESRI enterprise geodatabases on SQLServer: A and B. They are almost the same database (structure)
I have this code snippet in python to deleted several rows from versions:
def delete_project_features(edit_workspace, delete_tables, version, inspection_id):
    sde_conn = arcpy.ArcSDESQLExecute(edit_workspace)

    sde_conn.execute("EXEC sde.set_current_version '{0}';".format(version))
    sde_conn.execute("EXEC sde.edit_version '{0}',1 WITH RESULT SETS NONE;".format(version))

    for table in delete_tables:
        inspectionId = inspection_id
        feature_count = sde_conn.execute("SELECT COUNT(*) FROM {0}_evw WHERE INSPECTION_ID = '{1}';".format(table, inspectionId))
        sde_conn.execute("DELETE FROM {0}_evw WHERE INSPECTION_ID = '{1}';".format(table, inspectionId))
    try:
        sde_conn.execute("EXEC sde.edit_version '{0}',2 WITH RESULT SETS NONE;".format(version))

    except Exception as e:
        if repr(e).find("error 266") > -1:
            logEngine.log(_("Deleted all redundant features"), logger, SERVICE)
        else:
            mess = "Failed to delete features"
            raise Exception(mess)

When I run it on Database A, the rows are deleted as expected, but in Database B the rows remain there.
I always get the error:
AttributeError: ArcSDESQLExecute: StreamExecute ArcSDE Extended error 266 [Microsoft][ODBC Driver 17 for SQL Server][SQL Server]Transaction count after EXECUTE indicates a mismatching number of BEGIN and COMMIT statements. Previous count = 1, current count = 0.

despite the fact of the error, the rows are deleted on A
What setting could be different in the databases that prevents in B the successful result of the script execution?
Noted differences so far:
Database  A:
Compatibility level                 2016,
Recovery model                      Simple,
Collation                 SQL_Latin1_General_CP1_CI_AS.
Database B:
Compatibility level  2012,
Recovery model       Full,
Collation             Latin1_General_CI_AS.
EDIT:
The SQL commands run with any issue if they are run from SSMS.
The EGDB release for both is 10.7


texts: How to configure AGOL popup to open local files (PDF\JPEG)?

texts:
I'm trying to enrich my AGOL application by adding a hyperlink to ancillary documents that I have, on a permanent path at my local machine.

The question is this:

Is it possible to configure AGOL popup with HTML code to enable open files from local machine?

I've tried adding the path based on a the software's configuration wizard. Trying to insert an additional HTML prefix:

Click here to download

Doing so makes the browser open the path with my AGOL home page URL before the path and hence the document isn't opening.
It seems that ESRI blocked this option in order to increase credits use.


texts: Trying to Add Field and Calculate Area of a selected shapefile

texts:
I'm a Python beginner working with arcpy. I want to take a selected shapefile and add an area field to it, and then populate that field.

As of now, my code only works if I have the selected layer as a parameter input by the user. If I substitute that for the "selected layer" (which is commented out right now), the field is created, but it does not populate. I cannot figure out why.

What do I need to do in order to make this work?

import time, arcpy, sys, os, pythonaddins

## Set parameters

mxd = arcpy.mapping.MapDocument("CURRENT")
df = arcpy.mapping.ListDataFrames(mxd)[0]
mxdName = os.path.basename(mxd.filePath)
directory = os.path.dirname(mxd.filePath)
layers = arcpy.mapping.ListLayers(mxd, "*", df)
arcpy.env.overwriteOutput = False
#####SelectedLayer = pythonaddins.GetSelectedTOCLayerOrDataFrame()
SelectedLayer = arcpy.GetParameterAsText(0)
f_name = "Area_ha"
f_type = "DOUBLE"
f_precision = 5
f_scale = ""
f_length = ""
f_alias = ""
f_null = "NULLABLE"
f_required = "NON_REQUIRED"
f_domain = ""

## Create Area field, calculate area in hectares

arcpy.AddField_management(SelectedLayer, f_name, f_type, f_precision,
f_scale, f_length, f_alias, f_null, f_required, f_domain)
areaexpression1 = "{0}".format("!SHAPE.area@HECTARES!")
arcpy.CalculateField_management(SelectedLayer, f_name, areaexpression1,
"PYTHON_9.3", "")


texts: Converting from kilometers to degrees distances?

texts:
We have geometries in a Postgis database with SRID WGS84 and we have found lookups directly in degrees to be much faster than in kilometres, because the database can skip the projections I would think.

So, having a reference point in coordinates and a distance threshold in kilometres, is there a somewhat precise way to convert the distance to degrees for the lookup?

Example: Finding all locations within 10 Km of POINT(-3, 40). How do I convert the 10 Km to a degree value?


texts: Using ArcGIS Network Analyst to plot route for known OD pairs

texts:
I've recently been given a list of 'routes' from a client. The list simply comprised two sets of coordinates, for the start and end point of each route. In this dataset there are around 350 routes. I firstly plotted the start and end points as two separate feature classes.

My client wants to see the routes through the road network, rather than a simple crow-flies pairing. Using Network Analyst within ArcGIS 10.3, the best option seems to be to run a 'Closest Facility' calculation. I treat the start points of my routes as facilities, and the end points as incidents, and sort/differentiate both datasets using the route name field provided by the client.

As some of these routes overlap, and don't necessarily represent the closest pairing of points, I'm not actually looking for the closest facility at all, but this seems to be the only network analyst tool within the extension that produces actual routes through a network as an output. To get round this I've ramped up the number of 'Facilities to find' to 50.

Having run the calculation, the tool produces a series of polylines. I then have to export this as a new featureclass, and split the 'name' field (the result field). I then look for pairs where the first half of the 'name' field is equal to the second half, as the name field simply records the route names of the start and end points. However having run this tool I've found that some start points are never routed to their corresponding end points if they are particularly far away and (presumably) there are fifty closer end points, or incidents, for the software to analyse.

I suppose the obvious option is to cut my dataset down into a number of subsets. Having tried to increase the number of incidents for the tool to find my computer promptly ran out of memory (to paraphrase the error message).


texts: Working with lidar points (LAS To Multipoint) in ArcGIS for Desktop?

texts:
Currently I'm working with the LAS format and want to transform from WGS84 to RSO Kertau Malaysia.

I'm using ArcGIS-Las to Multipoint and Data Management Tool-Project and succeeded in the conversion.

However the output gave only 1 default class, as the raw files are already classified into their respective classes (ground, vegetation, buildings, etc.).

Are there any other ways of conversion which at the same time creates the classes?


texts: Setting transparency of layer group with Python in QGIS?

texts:
How can I set the transparency for all layers within a group of layers?

I assume this will require the use of a FOR loop in the Python Console.

Based on How to set the transparency for multiple layers or add a global transparancy preference?

I tried the following:

for layer in iface.legendInterface().layers():
   layer.renderer().setOpacity(0.5)

But I only get this error:

Traceback (most recent call last):
  File "<input>", line 2, in <module>
AttributeError: 'QgsVectorLayer' object has no attribute 'renderer'

What am I doing wrong?


texts: Two persons working on same layer in QGIS

texts: We are two people who need to edit on the same layer in QGIS. For example one of us will edit locally, the other one via a server connection.


texts: Copy tables from one server to another in PostgreSQL

texts: I am trying to copy a number of tables from one server to another using PostgreSQL as part of reorganising our data structure. I am pretty new to PostgreSQL so haven't much experience with it. Is there a simple way to do this?


texts: Dissolve with dominant/majority field using ArcGIS for Desktop?

texts:
I want to dissolve a feature class of land use polygons, joining subparcels in larger parcel polygons. I want the parcel polygon to keep the dominant land use of all subparcels, taking into account the surface that each subparcel takes.

Dissolve GP tool lets me calculate statistics over one field (min, max, avg) but I can't find a way to keep the dominant value.

I'm using ArcGIS 10.1


texts: LiDAR dataset displaying 500 meters to the east of correct position

texts:
Why does my LIDAR dataset show up ~500 meters to the east of the correct position even though a spatial reference/coordinate system is provided?

Some background:

The LIDAR data were downloaded from the USGS website as .zip files
Unzipped, the .las files have an associated .xml file with correct coordinates for tile position (UL Corner, UR Corner etc.)
Projected coordinate system (State Plane NAD 1983 Florida East Meters) is listed in .xml
The FL County layer in the image below is in the correct location.

How can I permanently bind the projection/spatial information to the .las file (not need the separate manually created .prj file)? I'm trying to use the LAStools (from lastools.org). When I use the LASclip, the clipping polygon is in the correct location but arcmap thinks the .las file is 500 meters east still, so nothing (or everything) gets clipped.


texts: Trouble working with a .TSV containing WKT

texts: I'm working with a .TSV file that represents district boundaries in a portion of a country in Carto. The polygon data is represented in a string column of Well-Known Text lat/long vertices, which I'm having trouble geocoding and visualizing. How would I go about representing the district boundaries as polygons on a map? Here is a screenshot of the data set for your reference.


texts: Download data on click via Leaflet webmap

texts:
I am working on a webmap which will allow users to click on a portion of the globe and download associated river data.  Here is my code so far.
<!DOCTYPE html>

    <title>GRWL test</title>
        <link rel="stylesheet" href="https://unpkg.com/leaflet@1.3.1/dist/leaflet.css"
   integrity="sha512-Rksm5RenBEKSKFjgI3a41vrjkw4EVPlJ3+OiI65vTjIdo9brlAacEuKOiQ5OFh7cOI1bkDwLqdLw3Zg0cRJAAQ=="
   crossorigin=""/>

This code outputs a grid system over the basemap as illustrated below, each grid is its own feature with an identifying attribute.  How would I allow the user to click on a tile of their choosing and download a shapefile whose name matches the identifying attribute of the tile.  All shapefiles will exist in a folder on my machine, eventually the shapefiles will be pushed to a server to allow functionality to non local users.


texts: I have created a RGB composite in SAGA GIS using the RGB composite tool. A one band raster is created, which can be visualised in SAGA using the option RGB coded values. When I try to import this raster in QGIS, the image is loaded in grayscale, because the RGB values are stored as a single number, for example "6512354". In SAGA this numbers are recognized, but not in QGIS. Is there any way I can visualise  the RGB coded values in QGIS?


texts: I'm planning buying a tablet for running Ubuntu. There are little info on experience with Ubuntu Touch and QGIS and does it work? What I want is to run it on the tablet, drawing polygons etc on screen, but which version works best? Anyone that have good experiences? And no Windows is not an option, maybe as dualboot.


texts: Working out percentage in QGIS field calculator

texts:
I have a project in QGIS where one of the things I need to do is work from a police data file with 9 months of crime information, linked to details of county population (age, employment, etc) to report the % of different crimes in each county (7 counties). When I follow the example screenshot given, word for word, to work out a %, I get an message which says: "An error occured while evaluating the calculation string.".

I'm not sure where i'm going wrong. Any guidence would be a lifesaver.


texts: Changing PostgreSQL database for QGIS project?

texts:
I have a QGIS project which uses PGSQL in the back end.  I am now setting up multiple instances of the database (dev, working and published) with a mechanism to migrate data between  the in a systematic manner.

Is there a straight forward way of changing the database name in a QGIS project so I can quickly test things after performing migrations?  i.e just change the database leaving schema, tables. etc the same.


texts:
I am calculating the path distance in ArcMap 10.1,

My inputs are point locations(source) from which I am calculating the path distance.
A cost raster that shows the impedance per cell(30m) and impedance in terms of minutes taken to cross each cell.

A DEM defining the elevation values at each cell location that is the surface raster and the same DEM to apply for slope correction ToblerAway

When I run the function without the DEM for slope correction I get realistic time of 0-80 minutes.

When I incorporate the DEM and use the ToblerAway within the path distance, my results become unrealistic and get 0- 0.36 minutes. The DEM is SRTM 30 m projected in UTM.
 What could be the problem?


texts: Setting grid in QGIS that automatically adjust its intervals?

texts:
I am setting a template for multiple maps in QGIS 3 Layout and i was wondering if there is a way of automatically adjusting the grid's interval as a function of the scale of the map?

I did a bit of research here:

Is there any other option to display DMS grid and graticule in QGIS Print Composer?
How to create a reference grid in QGIS print composer?

but they seem not to answer my question.

I want to make sure that no matter the scale of my map, this template adjusts the grid. Currently, one grid may work well for a determined scale

But the grid may get quite dense if i zoom out (see numbers in the top as well)

Or useless if i zoom in

I want to define the intervals in a way that they are automatically updated as a function of the scale of the map.

How do I do this?


texts:
I’m using a while statement to count down from x until it = y. Each loop decreases x by 1. This variable (among other variables)  is passed into a geoprocessing tool’s  expression variable as a string. The string exactly mimics a user specified expression. Alas the ArcPy .gp does not accept my identical variable as an expression.

How does one take a text string, and pass it off as a text expression?

####### Raster values to separate rasters
import arcpy
arcpy.CheckOutExtension("spatial")

######## Script arguments

FolderLocation = "C:\\Data\\Testing\\"
InRaster = "C:\\Data\\ basin00_dem"

TopVal = 236
BotVal = 127
Range = TopVal-BotVal
print "There are: "+ str(Range) +" values to iterate"
print " "
SingleValue = "s1_"

MegaCount = TopVal

while (MegaCount) > BotVal:
    MegaCount = MegaCount-1
    ## Reclass Single Values
    SingleValueName = SingleValue + str(MegaCount) +".tif"
    SingleValueNamePath = FolderLocation + SingleValueName
    RC1 = "\"" + str(BotVal) +" "
    RC2 = MegaCount -1
    NODATA = " NODATA;"
    RC3 = str(MegaCount) +" "
    RC4 = str(TopVal)+ " "
    Val = "1 " + "\""
    Qu = "\""
    Equ = str(RC1) +str(RC2) + NODATA + str(RC3) + str(RC4) + str(Val) +str(Qu)
    Range = Range -1
    print "Number: " + str(Range) + " Equ: " +str(Equ)
    arcpy.gp.Reclassify_sa(InRaster, "Value", Equ , SingleValueNamePath, "NODATA")

# WORKING .GP
# arcpy.gp.Reclassify_sa(InRaster, "Value", "127 233 NODATA;234 236 1", SingleValueNamePath, "NODATA")
# Example print of Equ: "127 129 NODATA;130 236 1"


texts: Is it possible to wrap date line with gdal.VectorTranslate()?

texts:
I use '-wrapdateline' option in ogr2ogr to split the polygons at the anti meridian. But in ogr2ogr Python port bWrapDateline is hardcoded to False.

Is there a way to pass 'wrapdateline' option to gdal.VectorTranslate() to split a single polygon into a multipolygon? I'm using GDAL 2.1.0


texts: Maximum ndvi for 3 months average in Google Earth Engine

texts:
I have a ndvi series which I am able to plot using below code:

var ndvi = l8.map(function(image) {
  return image.select().addBands(image.normalizedDifference(['B5', 'B4']));
});
var ndviChart = ui.Chart.image.series(ndvi, point, ee.Reducer.mean(), 500);
print(ndviChart);

My goal is to calculate the average ndvi for every 3 months window and then calculate the maximum ndvi over that 3 months window. How to access data in ui.Chart.series and is there any way by which I can do the above calculation on Google server ?

In javascript it can be done using

<script>
var lista = [22,4,5,6,11];
var maxAvg = 0.0;
for (var i = 0; i < lista.length - 2; i++) {
var avg = 0.0;
for (var j = 0; j < 3; j++) {
   avg = avg + lista[i+j];
}
if (maxAvg < (avg/3.0)) {maxAvg = (avg/3.0)};

}

document.write(maxAvg);
</script>

I want to replace lista with ndviChart series. Can anybody help me in this direction?


texts: Select by location is empty when done through arcpy but returns values when done manually

texts:
I have nearly 1000+ shapefiles that I need to (1) reproject, (2) covert to a raster, and (3) clip.  In order to make sure that there are values within the extent of concern, I put in a SelectByLocation_management logical check.  This is where I am having issues.  When I do the SelectByLocation, it comes up blank through the script.  However, when I do it by hand through the Selection>Select by Location toolbar, it selects the locations where there is an intersection.

Any ideas?

import arcpy, os
... arcpy.env.extent = "-1980501.4337461 2059002.55365216 -1940501.4337461 2090502.55365216"
... arcpy.env.workspace = r"C:\\Data_Prep\\2013_Shape"
... shapefiles = arcpy.ListFeatureClasses()
... projectOut = r"C:\\Data_Prep\\2013_Shape\\Project"
... rasterOut = r"C:\\Data_Prep\\2013_Shape\\Raster_Data"
... Extent_Polygon = 'Extent_Polygon'
... layer_files = r"C:\\Data_Prep\\2013_Shape\\Layers"
... error_list = []
... out_of_extent = []
... successful = []
... for shape in shapefiles:
...     print shape
...     # Process: Project
...     shape_name = os.path.join(projectOut, shape)
...     print shape_name
...     arcpy.Project_management(shape, shape_name, "PROJCS['USA_Contiguous_Albers_Equal_Area_Conic_USGS_version',GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]],PROJECTION['Albers'],PARAMETER['False_Easting',0.0],PARAMETER['False_Northing',0.0],PARAMETER['Central_Meridian',-96.0],PARAMETER['Standard_Parallel_1',29.5],PARAMETER['Standard_Parallel_2',45.5],PARAMETER['Latitude_Of_Origin',23.0],UNIT['Meter',1.0]]", "", "GEOGCS['GCS_North_American_1983',DATUM['D_North_American_1983',SPHEROID['GRS_1980',6378137.0,298.257222101]],PRIMEM['Greenwich',0.0],UNIT['Degree',0.0174532925199433]]", "NO_PRESERVE_SHAPE", "", "NO_VERTICAL")
...     print "%s was projected" %shape
...     #determining if any cells are within extent
...     layernm = shape[11:-4]
...     layer = arcpy.MakeFeatureLayer_management(shape,layernm)
...     arcpy.SelectLayerByLocation_management(layer, "INTERSECT", "Extent_Polygon","","NEW_SELECTION")
...     result = arcpy.GetCount_management(shape)
...     count = int(result.getOutput(0))
...     if count > 0:


texts: PyQGIS: Inserting a logo image in the layout

texts:
I have a big logo image (1443 x 1453) that needs to be imported and added to the print layout. The image should be resized (to about 15% of its original size) and placed at a specific location in the layout.
What would be the command for such an action?


texts: gdal_grid out of memory allocating 64 bytes

texts:
Python 2.7.4 32-bit
GDAL 1.10.1

I have a Python script which takes a point file with x,y,z coordinates, makes a vrt of it and then tries to create a Geotiff from it with gdal_grid.
During testing all goes well and then I try to put in a 3.5 GB file. It fails.
I try a junction of that file and it works but at a certain moment when my test file gets to big it fails again. The error is always the same:

ERROR 1: CPLMalloc(): Out of memory allocating 64 bytes.
FATAL: CPLMalloc(): Out of memory allocating 64 bytes.
ERROR 2: Cannot allocate 6563 bytes

The last amount of bytes is equal to the width of the geotiff it should create.
The fact that it works with smaller files makes me think I'm running into restrictions of the software.

If that is the case I would like to know if there is a workaround?


texts: arcpy: creating polyline features from the following text file format

texts:
I am looking into a way of reading the following file format with python, where

0 0, 0 1, 0 2, ect... is the polyline object number and 0(1) 464736.1645 260940. 098377 are x, y coordinates of polyline feature start/end.

How I can scan through this text file to read a polyline object and correspsonding start/end coordinates.

PS. CreateFeaturesFromTextFile ESRI script tool does not work properly in arcpy (arcgis 10) hence the need to create polylines from text file.

Your help is much appreciated.

Many thanks,

Meryloo


texts: Merge vector layers to fill gaps in QGIS

texts:
I have a zip code map which I added per zip code data too (a regional director for each zipcode), and then dissolved the map based on the director attribute.  As you may or may not know, zip codes are really just postal routes, so areas in the U.S. that do not receive post, do not have zipcodes.  Because of this, there are empty spaces on my map I want to fill.

So my question is, can I somehow merge my dissolved zip code map with a map of the US state boundaries, and use that to fill the empty spaces?  I guess I could manually draw out the polygons to fill the gaps, but this would be incredibly time consuming.

Here is a snapshot of Maine to show what I'm explaining.  Cyan is zipcode vector, grey is US state boundaries.

Any ideas?


texts: Hosting WMS layer using Geoserver, Tomcat, Amazon AWS?

texts:
My goal: provide clients with a password protected WMS url so they can access some of our GIS data.

Where I'm at so far:

OS: windows 8

I have a polygon shapefile created in QGIS I want to share as a WMS. This shapfile has been uploaded to a PostgreSQL database with the PostGIS extension enabled on the schema. I use PGadmin 4 to do this. The PostgreSQL database is stored on Amazon AWS.

Our client who is to receive the WMS layer, has been added as a read-only user role type. I can successfully log in to the database using the credentials I created for them.

I have installed tomcat apache with geoserver installed on it (127.x.x.x:82/geoserver). In geoserver (via tomcat web application manager) I create the WMS layer and style it. After publishing the layer as a WMS I can view it in QGIS.

However, on another PC, I can't access the WMS layer at all. So it seems the WMS is still hosted locally on the original PC. Even though I have tomcat installed and I'm not using localhost.

Now, I'm a bit of a noob with web hosting etc, so what should I do to gain access to this WMS layer on other PCs (online instead of local)?


texts:
I have part of a topographic map with contour lines and height points, but I don't have any further info about it. I know it must be in Balkan Zone 7 CRS and nothing more.

How should I georeference it (in QGIS)?


texts: How to draw a Polygon in QGIS 3.8 on Mac?

texts:
My apologies, this is quite certainly a duplicate of an old question, but I cannot find the answer. I am struggling with drawing polygons in the new QGIS 3.8.3 on Mac.

In QGIS 2 drawing polygons was as simple as enabling the editing mode, clicking the "Add Polygon" button, drawing the polygon by clicking and then right-clicking (Ctrl + click on the trackpad). After this, a pop-up would come and ask to add Feature Attributes.

In QGIS 3.8 everything works relatively similar until the right-click. Nothing happens. Only Esc works and removes the polygon (as it should). Which button am I supposed to press?


texts: Polygon to Raster conversion: Raster shifted ~10m horizontally. How to fix?

texts:
I converted a number of polygons into a raster. The cell size is correct but the raster has shifted ~10m horizontally.

Is there an issue with my method or is there a tool I can use to fix this?


texts: What's the difference between NAD83(2011) (OGR WKID 6318) and GCS_NAD_1983_2011 (ESRI WKID 6318)?

texts:
I just projected a feature class originally in Web Mercator (WKID 3857) to WKID 6318 (ESRI calls this "GCS_NAD_1983_2011") using ESRI's projection tool, and the same feature class to WKID 6318 (OGR calls this "NAD83(2011)") using ogr2ogr. They are very different:  if I put them side by side in ArcMap they are 38" off in some places.  What's more, I added the GCS_NAD_1983_2011 shapefile (the ESRI-projected shapefile) to ArcMap first and when adding the NAD83(2011) shapefile (the ogr-projected shapefile) ESRI prompted to use a transformation since they are different datums.

This is unexpected:  I expected WKID 6318 to be the same between the two (ESRI and ogr).  Why would they be so different?

This is the full description of GCS_NAD_1983_2011/ESRI (looking at the 'XY Coordinate System' tab in the Properties window of ArcCatalog):
GCS_NAD_1983_2011
WKID: 6318 Authority: EPSG

Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_NAD_1983_2011
  Spheroid: GRS_1980
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314140356
    Inverse Flattening: 298.257222101

And this is the same info for NAD83(2011)/ogr:
NAD83(2011)
Authority: Custom

Angular Unit: Degree (0.0174532925199433)
Prime Meridian: Greenwich (0.0)
Datum: D_NAD83_National_Spatial_Reference_System_2011
  Spheroid: GRS_1980
    Semimajor Axis: 6378137.0
    Semiminor Axis: 6356752.314140356
    Inverse Flattening: 298.257222101

Everything is the same except for 'Datum'.

Why would that be?  Which one is the "authoritative" datum??


texts: Export GRASS raster as GeoTIFF : nodata values exported as -nan (what happens in R?)

texts:
I am using the package raster to work with GeoTIFF raster files which I exported from GRASS using r.out.gdal. GRASS assigned the value -nan to some raster files (type FLOAT32) for missing values. Other files had 65535 or 255 assigned to missing values.

In R, I reassign these values to missing like this:

raster[raster==255] <- NA

The problem is that I do not know what to do with multiple rasters that had the -nan assigned to missing values. For what I see when I examine the structure of these rasters (which I have as part of a list), the -nan are stored as NaN.I tried doing:

rasterlist[["element_i"]][rasterlist[["element_i"]] == NaN] <- NA

but I do not see any difference, for example when I do:

getValues(rasterlist[["element_i"]], row=1000)

The output is something like:

8100 NaN 8101 NaN NaN 8102 8104 8106

so the NaN are not transformed into NA.

Any suggestion?


texts: Where do I find official summit elevations found on USGS topo maps?

texts:
I'm creating topo maps for my day hikes guides for national parks. I've already done one for Glacier National Park. When I was working on summits, I downloaded the data for Montana from GNIS and when I used the supplied elevations, I found it was different than the National Geographic maps, USGS topos and official Park handouts (all three had the same elevations). The problem appears to be that GNIS has DEM derived elevations.

I've scoured USGS's national map viewer, and the web at large trying to locate the data and am coming up short. I find pockets of data that seem to be grass roots efforts to basically do what I did and pull the elevations off of a USGS topo and put it in a table. I'd rather have the "official" elevations to save time and help with accuracy (instead of the propensity for me to fat finger it).

Anyone have a place to look?


texts: Splitting line with equal parts in range

texts: How to split line into equal parts in QGIS, with a certain condition: the parts to be between 50m to 150m?


texts: Adding Esri Elevation to QGIS, possible?

texts:
I can add World Imagery through Layer > Add Layer > Add ArcGIS MapServer Layer, using this URL. Is this somehow possible with Esri's World Elevation layer (info see here) as well?

When I try to add it as layer using this URL, QGIS does return a layer ("Elevation(90m/1km)"), but can't add it (error "Layer is not valid"). I guess it's because it's a GlobeServer?


texts: Combine multiple partially overlapping rasters into a single raster in R

texts:
I have about 1,000 rasters which are around 10km x 10km and partially overlap to cover a whole country.

I would like to convert these into a single national raster, and where raster overlap, add the cell values together.

Currently, I am using ArcGIS's workspace to new raster tool, but I would like to replicate the process in the R language.

I've tried a couple of things but usually end up with only the intersection of the rasters. Many of the solutions on StackExchange advocate cropping to a common area, which is the exact opposite of what I want.

IS there a general solution for aligning two rasters and then making a mosaic which has the combined extent of both rasters?

Can this method be scaled to work for thousands of rasters?

EDIT

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

locale:
[1] LC_COLLATE=English_United Kingdom.1252  LC_CTYPE=English_United Kingdom.1252    LC_MONETARY=English_United Kingdom.1252
[4] LC_NUMERIC=C                            LC_TIME=English_United Kingdom.1252

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] rgdal_1.2-5  raster_2.5-8 sp_1.2-4

loaded via a namespace (and not attached):
[1] tools_3.3.2     Rcpp_0.12.9     grid_3.3.2      lattice_0.20-34

> rasterOptions()
format        : raster
datatype      : FLT8S
overwrite     : FALSE
progress      : none
timer         : FALSE
chunksize     : 1e+07
maxmemory     : 1e+08
tmpdir        : D:/RTemp\RtmpgzISfL/raster/
tmptime       : 168
setfileext    : TRUE
tolerance     : 0.5
standardnames : TRUE
warn depracat.: TRUE
header        : none


texts: How to set style for a Wind - point layer on GeoServer

texts:
I am serving a Wind as a point layer on GeoServer.
In each point I have data:

wind speed [m/s] and
wind direction [degrees]

Currently I am using default style, but I would like to be able to change its representation to

Arrows (http://geoserver.geo-solutions.it/edu/en/multidim/accessing_multidim/rtx/currents.html) and/or
Wind Barbs (http://geoserver.geo-solutions.it/edu/en/multidim/accessing_multidim/rtx/wind_barbs.html), so I can use whatever visualization will be necessary.

Is it possible to do it with my data source and how can it be done?


texts: QGIS plugin builder puts icon in different toolbars depending on plugin template selected

texts:
When building a plugin using the QGIS Plugin Builder plugin, if I select 'Tool button with dialog' as the template, then my plugin button will be placed in the Plugins toolbar in QGIS.

If I build the plugin boilerplate with the 'Tool button with dock widget' template, then the plugin icon is put in its own toolbar.

I have compared the code and cannot find how to change this. I would like to build a plugin with a dialog, but have the icon in its own toolbar.

Please let me know if there is a way to make my question more clear.


texts:
Suppose the point has coordinates for a search within the box:
(long:-155, lat:19)

If we have the box: SW - (-170, 10), NE (29,69) that everything is fine and we can find this point because it's inside this box.

But if we move the box a little bit the coordinates are SW - (144, 10), NE (-75,69)   that results are zero although the point is still inside of this box.

How to do it properly to find the point in the second situation?


texts: OSS java approach to convert Shapefile to GML

texts:
We've got a project coming up where we need to upload  .zip of a shapefile, join the .dbf to some lookup tables, join to other tables also uploaded (one to many attributes for each feature), and ultimately generate GML following a defined XML schema. So the geometry in the GML will come from the Shapefile (lines and polygons) and the attribs will come from the core table and some additional uploaded tables.

No projection needed, no searching, no filtering. The system is a Java web app and we need to give the source out as open source when we are done (no proprietary code). We just need something simple, easy and modern.

Would you use Geotools for this or call ogr2ogr? Is geotools overkill in terms of complexity for this?


texts: How do I declare a QGIS custom expression function and use it in a single python script?

texts:
I would like to create and populate a virtual field in a QGIS vector layer attribute table using a python script.  I think the easiest way to do this is to create a custom expression function, and, then use this function to populate the virtual field.  How do I instantiate a custom expression function from within a python script?  What I have tried with Python 3 from a Python IDE is:

    from qgis.core import *

@qgsfunction(args=0, group='Custom', usesGeometry=False)
def a_function():
    return "'Test'"

QgsExpression.registerFunction(a_function)
print(QgsExpression(a_function()).evaluate())

This returns a TypeError:

print(QgsExpression(a_function()).evaluate())
TypeError: 'QgsPyExpressionFunction' object is not callable


texts: Get the current serial protocol config of a TU30-D400 Series GPS

texts:
I have an old Jupiter GPS from the TU30-D400 Series

When I put its 7th to the GND, I can force it to use the NMEA Protocol with a 4800 baud rate.

The thing is that GPS is used by a pretty old tide gauge, and it has its own serial configuration.

How can I know what kind of config does it have when I'm on the NMEA Protocol ( Pin 7 to GND )

http://www.sander-electronic.de/datasheet/gps/NAVMAN_Jupiter_11_Datasheet.pdf


texts: Looking for the best data source for a geocoding database for Canada

texts:
I'm planning to build a Geocoding database for Canada. The goal is that - the database should be able to provide latitude and longitude for any given street level address in Canada. A similar implementations for USA (Geocoder::US 2.0, jGeocoder ) uses US Census Data (Tiger).

Provincial Base Mapping Data Sources in Canada provides a list of data sources for Canada (both nation wide and provincial based). I couldn't find anywhere is a discussion about the quality of these sources. Also I don't have much time to go through all of them in detail and pick the best.

If anyone had done similar research before about various geo-spatial/geo location sources of Canada for similar implementation. Or if found any similar implementations/attempts some where on web. Guide me please. Any suggestions regarding my goal are also appreciated.

Note: I don't want any suggestions for services that can help like Google API, yahoo place finder etc.


texts: Updating a table with the changes made from another table

texts:
I am trying to do updates to a source table from another table that has a like field (propertyid). I am trying to update all of the fields in the source table. I have gotten the outcome I want with ModelBuilder, but the model is inefficient because I would have to add a Field Calculator tool for every field that I am trying to update. Here is what I have written, it is running but nothing happens in my output.

import arcpy

# Define workspace

arcpy.env.workspace = "C:\Users\CartoFront\Desktop\Updates\Updates.gdb"

#Define Feature Classes
fc1 = "Update_121316_test"
fc2 = "ChiMetroJoin12012016clean"

#Define Fields
fc1_fields = ["*", "propertyid"]
fc2_fields = ["*","propertyid"]

#Define Values
fc1_value= "*"
fc2_value = "*"

cur = arcpy.da.SearchCursor(fc1, fc1_fields)

for row in cur:
    cur2 = arcpy.da.UpdateCursor(fc2, fc2_fields)
    for row2 in cur2:
        if row[1] == row2[1]:
            if row[0] == fc1_value:
                row2[0] = fc2_value
                cur2.updateRow(row2)
del row, cur, row2, cur2