texts: Converting .tif files to vector MBTiles by using gdal_polygonize

texts:
I am trying to use GDAL 2.2.2 and the gdal_polygonizecommand to convert a TIF into vector MBtiles from the link:
https://github.com/geostarters/gdal_snipets
I run the following command on the tif as in the link:
gdal_polygonize.py -8 -nomask -b 1 input_4326.tif \
-f "MBTiles" output.mbtiles \
-lco NAME=nameoutput \
-dsco WRITE_BOUNDS=YES \
-dsco MAX_FEATURES=20000000 \
-dsco MAXZOOM=15

However I get the following ERROR:

Traceback (most recent call last):   File
"/usr/bin/gdal_polygonize.py", line 180, in
dst_ds = drv.CreateDataSource( dst_filename )
AttributeError: 'NoneType' object has no attribute 'CreateDataSource'

It seems to me that I do not have the GDAL MBTiles Driver 'installed':
https://gdal.org/drivers/raster/mbtiles.html
I installed gdal onto Uubntu 16.04 via the terminal.
Is there anything I can do to get this working?


texts:
QGIS Zurich 3.18.2 -Win10
I have a virtual layer that has multiple polygons each with two pertinent Attribute fields. "Field1Name" and "Field2Name". I can set and Style the Label for one field and I can display these two fields as a Label, using the Expression Value "field1name" || '\n' || "field2name".
BUT is it possible to set a different Text/Render Style for each of the two fieldnames in the Label to display at the same time? I want both fields data to appear within the polygon label, but one highlighted, bold and large and the other subordinated, smaller less obvious. If I set up two Text Styles up is there an Expression that will achieve this two component styled Label?


texts: What is the Pro or Cons for mosaic 30 GB jp2000?

texts:
Has anyone tried mosaic jp2000 imagery in ERDAS, ArcGIS , ENVI or any softwares that they have accomplished ?

I have three jp2000 files that is total about 28.1 GB that I want to mosaic with ?


texts:
I bet this is probably an issue with my PC, but I have a really odd error with ArcGIS lately.

Refer to the image, it'll explain it best:

Notice how their no map displayed? No matter what I do I can't get it to display. If I add data, change the basemap, restart my computer, re-installed Explorer. This started happening about a week ago for no apparent reason. Nothing had crashed or changed recently on my PC.

Has anyone ever heard of any bugs that might be the cause of this?

My ArcGIS Explorer build is: 1750 (2.0.0.1750)

Any help would be greatly appreciated.


texts: Is it possible to include the Map composer in Sextante for QGIS 2.0?

texts:
I would like to solve a problem (Species atlas (how to automate to create on page for each species with QGIS)).
To solve it by an other way, I would like to know if it is possible to make a sextante process, and at the end, call, or use the map composer, directly into the Models Builder?
If yes, how to do it?

Thanks.


texts: Splitting line vector into components in Google Earth Engine

texts:
I’m trying to split a line vector (with buffer, although that could come after I split it) into many different components each 500 meters long to use for analysis.
Is there an easy way to do this so I dont have to do it by hand?


texts: Extracting value of Landsat Band 8 at given location in Google Earth Engine

texts:
I would like to have the value of the B8 of Landsat at a given location pt as an output but I keep on getting errors. I think I confused some values but cannot find the mistake. My code is below.
// LANDSAT

// This example demonstrates the use of the Landsat 8 QA band to mask clouds.

// Function to mask clouds using the quality band of Landsat 8.var maskL8 = function(image) {
  var qa = image.select('BQA');
  /// Check that the cloud bit is off.
  // See https://www.usgs.gov/land-resources/nli/landsat/landsat-collection-1-level-1-quality-assessment-band
  var mask = qa.bitwiseAnd(1 << 4).eq(0);
  return image.updateMask(mask);
}

var start = '2019-06-01'
var finish = '2019-06-30'
var pt = ee.Geometry.Point([-49.31582, 69.56833]);
// Map the function over one year of Landsat 8 TOA data and take the median.
var l8 = ee.ImageCollection('LANDSAT/LC08/C01/T1_TOA')
.filterDate (start,finish)
.filterBounds(pt)
.map(maskL8)
.select ('B8');

Map.centerObject(pt, 16);

// B8 VALUE OF PT

var myB8 = l8.select("B8");
print("myB8",myB8);

var getB8 = function(image) {

  // Reducing region and getting value
  var value_B8 = ee.Image(image)
    .reduceRegion(ee.Reducer.first(), pt)
    .get('B8');

  return value_B8;
};

var count = myB8.size();

var B8_list = myB8.toList(count).map(getB8);

print("B8 list", B8_list);

var img = myB8.first();

Map.addLayer(img);
Map.addLayer(pt);

var allDates = myB8.aggregate_array('system:time_start');

var allDatesSimple = allDates.map(function procDates (ele) {

  return ee.Date(ele).format().slice(0,10);

});

var paired = allDatesSimple.zip(B8_list);

print (paired);


texts: ArcGIS Network Analyst not recognizing One-way Restrictions?

texts:
I have created a network and have been attributing my various routes in a domain with a one-way field. At current all values are From-To (FT) in the one-way field. However the network does not honor the one-way restrictions.

I figure its within the VB script but I cannot work it out.

For example the road below has a One-Way set for FT  but is not honored nor does the network recognize it if using Network Identify Tool.

Network identify tool showing transferable in both directions

I am also using the default values Arc provides upon recognizing a field, as seen below:

From-To

To-From

Both source data and network edges are running in the correct direction (shown by arrow at end symbology). Oneway attribute is set to FT, but the network still shows its completely traversable.


texts:
I'm using QGIS 2.6.1 and i've been having problems using the polygonize function.

I keep getting a "DLL load failed" import error. I've searched older threads which unfortunately weren't able to help. I tried resetting the GDAL executable path to C:/OSGeo4W/bin/osgPlugins-3.0.0, and the python path to: C:/Program Files/QGIS Brighton/apps/Python27/Lib/site-packages/osgeo.

I don't get this problem when I try to use other functions (Rasterize or Extraction-Contour, for example). So i'm not sure whether it's a problem with the path or something else. Any ideas?

  Traceback (most recent call last):
  File "C:\Program Files\QGISBR~1\bin\gdal_polygonize.py", line 36, in
  import gdal, ogr, osr
  File "C:\Program Files\QGIS Brighton\apps\Python27\Lib\site-packages\osgeo\gdal.py", line 26, in
  _gdal = swig_import_helper()
  File "C:\Program Files\QGIS Brighton\apps\Python27\Lib\site-packages\osgeo\gdal.py", line 22, in swig_import_helper
  _mod = imp.load_module('_gdal', fp, pathname, description)
  ImportError: DLL load failed: Le module spécifié est introuvable.


texts:
I have a layer of AIS points that I want to make into lines using st_makeline and then use st_intersects to determine if any of the lines intersect a polygon.

So I think I want to nest

st_intersects(st_makeline, polgon.geom)

The SQL I am using for the makeline is:

SELECT mmsi, ST_MakeLine(geom ORDER BY date_time_stamp) As newgeom
    FROM ais.cargo As ais
    GROUP BY mmsi;

How do I wrap the st_intersects around this?


texts: ArcGIS Server 10.5.1 timeouts/500 errors

texts:
We have 3 app servers running ArcGIS server 10.5.1, all three servers connect to a file server. The config store is on the file server, and that contains all the server and service content.

Intermittently, one of the servers will start error-ing out, not finding the service content, but oddly the ArcGIS Server service will be running on the problem server as well all the ArcSOC.exes.

We have to manually restart the ArcGIS server windows service, and then problem is resolved.

Has anyone encountered these symptoms before? What was the cause? What was the solution? I may have to just uninstall and reinstall ArcGIS server, but the odd thing is this installation has been running since January 2018 without issues until the last few months.


texts:
I've only just started using drawing tools in ArcGIS, and have the following question I haven't found the answer to online.

Is there a way of taking a drawn ellipse (at a slight angle) and adding say 50-100 of the same ellipse, arranging them as efficiently as possible so that there's as little blank space as possible around them (but no overlaps)?

I have tried doing this by hand but I can't get the snapping tools to work properly (plus it's pretty time consuming), so ideally looking for a way to get this done quickly and introducing little/no human error into the process.


texts: ValueError: Invalid field type  when exporting shp?

texts:
I try to export an shp from Python using geopandas with this:

#Create a output path for the data
out = r"C:\Users\hhh\Desktop\ge.shp"

## Write those rows into a new Shapefile
df.to_file(out)

and I get this error:

ValueError: Invalid field type <class 'shapely.geometry.polygon.Polygon'>

The type of the df is geopandas.geodataframe.GeoDataFrame and I have done some geoprocessing which may be causing errors at the output. Also, it gives an output but when I open it, it shows nothing and the attribute table is also empty.

I checked this with another case by taking a subset from an shp read in python:

out = r"C:\Users\hhh\Desktop\hid.shp"

selection = data[0:50]

selection.to_file(out)

and it worked.

What can be causing this problem in the first case?


texts: Understanding Projection for Near tool in ArcMap?

texts:
I need to find the nearest distance of points from points in ArcMap. Both layers are in WGS 1984.

When I choose geodesic distance in Near tool, I get results in meters.

I don't do any projection.

Is this wrong?


texts: PostGIS function to find zip codes with radius in miles?

texts:
I have very little knowledge and understanding Postgres/PostGIS. I am having trouble by using the ST_PointInsideCircle function for radius that uses miles or meters or degrees.
I think this function expects degrees in radius parameter. I have a table with all zip codes with lat/long and I have input with lat/long/radius (miles).
I need to select all the ZIP codes found inside that radius and provided lat/long.

Can anyone help me to find right function or script?

Postgres 9.5 /PostGIS installed.

lat=37.6210700000000031,
long=-122.429737000000003,
radius=5 miles

my code :

SELECT zip_code from mytable
  where ST_PointInsideCircle(ST_Point(mytable.latitude,mytable.longitude), 37.6210700000000031, -122.429737000000003, 5);


texts:
I set up a layer with "time" dimension in GeoServer, but I'm having quite a hard time figuring out how to use it...

It is a vector layer with some polygons in it. This is my layer configuration for Time Dimension:

Since I set a valid default value, I am able to get that default value simply calling my layer:

http:///geoserver/emodnet/wms?service=WMS&version=1.1.0&request=GetMap&layers=emodnet:SatData&styles=&bbox=-180.0,-90.0,180.0,90.0&width=660&height=330&srs=EPSG:4326&format=application/openlayers&

but I can't possibly figure out how to get the map for different time values!

Also, is it possible to use the 'time' parameter in the WFS too?


texts: How to disable QGIS asking for DB password

texts: I open a QGIS project and I am prompted for underlying database password. For certain reasons I do not want to connect to underlying database (let's say I just want to change a certain stored local shape layer style) so I click "Cancel". However almost on every movement in the desktop application I am prompted with this dialog. How to turn it off for good?


texts: How can I programmatically get the path of "Python.exe" used by ArcMap

texts:
I am working with an add-in of ArcMap in C#. From C# code, I have executed some Python scripts. Now, to run those script, I have hard-coded python path. But this is not portable. So, I want to get the path of the Python executable from the code and use it.

Question:

How can I get the path of the Python executable used by ArcMap from C# code?

EDIT :

From your suggestions, for now I am using "path environment" to get Python path.

//get python path from environtment variable
string GetPythonPath()
{
    IDictionary environmentVariables = Environment.GetEnvironmentVariables();
    string pathVariable = environmentVariables["Path"] as string;
    if (pathVariable != null)
    {
        string[] allPaths = pathVariable.Split(';');
        foreach (var path in allPaths)
        {
            string pythonPathFromEnv = path + "\\python.exe";
            if (File.Exists(pythonPathFromEnv))
                return pythonPathFromEnv;
        }
    }
}

But there is a problem:

When different version of python is installed in my machine, there is no guarantee that, the "python.exe" I am using, ArcGIS also using that.

I don't appreciate using another tool to get "python.exe" path. So, I really think if there any way to get the path from registry key. For "ArcGIS10.0" registry looks like:

And for that, I am thinking about following way to get the path:

//get python path from registry key
string GetPythonPath()
{
    const string regKey = "Python";
    string pythonPath = null;
    try
    {
        RegistryKey registryKey = Registry.LocalMachine;
        RegistryKey subKey = registryKey.OpenSubKey("SOFTWARE");
        if (subKey == null)
            return null;

        RegistryKey esriKey = subKey.OpenSubKey("ESRI");
        if (esriKey == null)
            return null;

        string[] subkeyNames = esriKey.GetSubKeyNames();//get all keys under "ESRI" key
        int index = -1;
     /*"Python" key contains arcgis version no in its name. So, the key name may be
     varied version to version. For ArcGIS10.0, key name is: "Python10.0". So, from
     here I can get ArcGIS version also*/
        for (int i = 0; i < subkeyNames.Length; i++)
        {
            if (subkeyNames[i].Contains("Python"))
            {
                index = i;
                break;
            }
        }
        if(index < 0)
            return null;
        RegistryKey pythonKey = esriKey.OpenSubKey(subkeyNames[index]);

        string arcgisVersion = subkeyNames[index].Remove(0, 6); //remove "python" and get the version
        var pythonValue = pythonKey.GetValue("Python") as string;
        if (pythonValue != "True")//I guessed the true value for python says python is installed with ArcGIS.
            return;

        var pythonDirectory = pythonKey.GetValue("PythonDir") as string;
        if (pythonDirectory != null && Directory.Exists(pythonDirectory))
        {
            string pythonPathFromReg = pythonDirectory + "ArcGIS" + arcgisVersion + "\\python.exe";
            if (File.Exists(pythonPathFromReg))
                pythonPath = pythonPathFromReg;
        }
    }
    catch (Exception e)
    {
        MessageBox.Show(e + "\r\nReading registry " + regKey.ToUpper());
        pythonPath = null;
    }
    return pythonPath ;
}

But before using the second procedure, I need to be sure about my guesses. Guesses are:

the "True" associated with python means python is installed with ArcGIS
ArcGIS 10.0 and upper version's registry key will be written in same process.

Please help me to get any clarification about my guesses.


texts: Projecting ArcGIS Tiled Service for Esri-Leaflet?

texts:
I'm following this example from the esri-leaflet docs.

I have an ArcGIS Tiled MapService in ESRI:102696. So I'm trying to use Proj4Leaflet to project the service and use it in my application.

I found some proj4defs here and here. They do differ slightly, but neither seems to work.

This is the Map Service I am trying to use.

I put together a JSBIN example with the code below.

  // The two proj4defs I found
  //var crs = new L.Proj.CRS('ESRI:102696', '+proj=tmerc +lat_0=35.83333333333334 +lon_0=-90.5 +k=0.999933 +x_0=250000 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192  no_defs', {
  var crs = new L.Proj.CRS('ESRI:102696', '+proj=tmerc +lat_0=35.83333333333334 +lon_0=-90.5 +k=0.9999333333333333 +x_0=250000 +y_0=0 +datum=NAD83 +units=us-ft +no_defs', {
    transformation: new L.Transformation(1, 17627100, -1, 19798600),
    resolutions: [
      303.8194444444444,
      217.01388888888889,
      86.80555555555556,
      65.10416666666666,
      43.40277777777778,
      20.833333333333332,
      8.680555555555555,
      4.166666666666666,
      2.083333333333333,
      1.0416666666666665,
      0.5208333333333333
    ]
  });
  // pass crs into map constructor
  var map = new L.Map('map', {
    crs: crs
  }).setView([38.5432816667, -90.2688116667], 5);  L.esri.tiledMapLayer('http://maps.stlouisco.com/ArcGIS/rest/services/Maps/PublicBase/MapServer', {
    maxZoom: 10,
    minZoom: 0,
    continuousWorld: true
  }).addTo(map);

I'm not sure if the transformation is off or what.


texts: crop wrld_simpl gives error after R update

texts:
This code was working for me before I updated R to 4.0, but now it throws an error. I get this error on a Mac and a PC. I updated gdal on my mac, but this did not help. I made sure the rgeos package is up to date. It does not seem to be related to the extent I choose (i.e. extent(-180, 180, 50, 90) throws the same error).
data("wrld_simpl", package = "maptools")
out <- crop(wrld_simpl, extent(130, 180, 40, 70))

x[i, ] is invalid
Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, unaryUnion_if_byid_false,  :
  TopologyException: Input geom 0 is invalid: Ring Self-intersection at or near point 143.66192817999999 49.312211990000002 at 143.66192817999999 49.312211990000002
In addition: Warning messages:
1: In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Ring Self-intersection at or near point 143.66192817999999 49.312211990000002
2: In rgeos::gIntersection(x[i, ], y, byid = TRUE, drop_lower_td = TRUE) :
  Invalid objects found; consider using set_RGEOS_CheckValidity(2L)


texts: Looking for XSL templates using Google Charts API?

texts:
I am searching for examples on usage of ArcGIS Desktop HTML popups in general, being more specific on xsl templates using google charts api on it.

Are there any references available?


texts: Intersecting geometry sum - sharing the boundary count as intersection

texts:
This question stems from Calculating proportion of area covered by another layer - not just the intersection

I have narrowed down the problem.
I am attempting to find the percent of district layer(layer D) covered by another layer B.
In the picture layer B consist of 1+2+3.
When creating the intersection between layer D and B, it splits them into 1,2,3.
I create the interaction, find area for each interaction, and in the district layer, create
100* intersecting_geom_sum('Intersection','area_intersection')/$area
The ideal end result is area of 1 / area of D.
But it calculates area of 1+2+3/ are of D.
I have checked that the geometries are valid, but I think the interaction command is the problem.
When sharing the boundary counts as intersection.
How do I restrict it such that sharing boundary does not count as intersection?


texts: PostGIS adding point to an existing multipoint

texts:
I use PostGIS and I would like to add a new point to an existing MultiPoint using ST_Collect.
The problem is that I got this error message:

  ERROR:  Geometry type (GeometryCollection) does not match column type
  (MultiPoint) CONTEXT:  PL/pgSQL function inline_code_block line 10 at
  assignment SQL state: 22023

What do you think? What is wrong in my code? The first ST_Collect runs successfully, but the second one is not.

My code:

DO $$
<<first_block>>
DECLARE
  geom Geometry (MultiPoint,0);
  astext varchar;
BEGIN
geom=ST_Collect(geom,ST_Multi(ST_GeomFromText('POINT(1 4)')));
geom=ST_Collect(geom,ST_Multi(ST_GeomFromText('POINT(1 10)')));
SELECT ST_AsText(geom) INTO astext;
RAISE NOTICE '%',astext;
END first_block $$;

I tried with GeometryCollection instead of MultiPoint:

DO $$
<<first_block>>
DECLARE
  geom Geometry (GeometryCollection,0);
  astext varchar;
BEGIN
geom=ST_Collect(geom,ST_Multi(ST_GeomFromText('POINT(1 4)')));
geom=ST_Collect(geom,ST_Multi(ST_GeomFromText('POINT(1 10)')));
SELECT ST_AsText(geom) INTO astext;
RAISE NOTICE '%',astext;
END first_block $$;

And I've got this error:

  ERROR:  Geometry type (MultiPoint) does not match column type
  (GeometryCollection) CONTEXT:  PL/pgSQL function inline_code_block
  line 9 at assignment SQL state: 22023

I tried with Multipoints too:

DO $$
<<first_block>>
DECLARE
  geom Geometry (MultiPoint,0);
  astext varchar;
BEGIN
geom=ST_Collect(geom,(ST_GeomFromText('MULTIPOINT(1 4)')));
geom=ST_Collect(geom,(ST_GeomFromText('MULTIPOINT(1 10)')));
SELECT ST_AsText(geom) INTO astext;
RAISE NOTICE '%',astext;
END first_block $$;

Error message:

  ERROR:  Geometry type (GeometryCollection) does not match column type
  (MultiPoint) CONTEXT:  PL/pgSQL function inline_code_block line 10 at
  assignment SQL state: 22023


texts:
i am using the c++ api qgis 2.8.2
i have a strange problem when i want to use the pantool for my canvas

i used this code :

mpPanTool = new QgsMapToolPan((QgsMapCanvas*)mpMapCanvas);
((QgsMapToolPan*)mpPanTool)->activate();
((QgsMapCanvas*)mpMapCanvas)->setMapTool((QgsMapToolPan*)mpPanTool);

building is good,
when i launch my application and i load a raster or a s57 map
if i click on the map or on the canvas nothing happen till i move the mouse
when i move the mouse it seems that the map change the upper left origin point and set a new one at the mouse position
by example if i left click on the map and i move the mouse, the map is automatically centered with the upper left point of the canvas at the mouse coordinate

perhaps i forget something, pan works but the problem is this coordinate change

please help


texts: Error trying to cast an ee.Image to a numpy array

texts:
I want to cast my ee.Image, to do so I followed this link  implementation, but I got this error:

Traceback (most recent call last):
  File "/home/felipe/.local/lib/python3.6/site-packages/ee/data.py", line 338, in _execute_cloud_call
    return call.execute(num_retries=num_retries)
  File "/home/felipe/.local/lib/python3.6/site-packages/googleapiclient/_helpers.py", line 134, in positional_wrapper
    return wrapped(*args, **kwargs)
  File "/home/felipe/.local/lib/python3.6/site-packages/googleapiclient/http.py", line 898, in execute
    raise HttpError(resp, content, uri=self.uri)
googleapiclient.errors.HttpError: <HttpError 400 when requesting https://earthengine.googleapis.com/v1alpha/projects/earthengine-legacy/value:compute?prettyPrint=false&alt=json returned "Image.sampleRectangle: Fully masked pixels / pixels outside of the image footprint when sampling band 'soilDiff' with no default value set. Note that calling sampleRectangle() on an image after ee.Image.clip() may result in a sampling bounding box outside the geometry passed to clip().">

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "exportImgs_v2.py", line 205, in <module>
    array = imageToArray(imgAlerta, rect)
  File "exportImgs_v2.py", line 41, in imageToArray
    print(band_arrs.getInfo())
  File "/home/felipe/.local/lib/python3.6/site-packages/ee/computedobject.py", line 95, in getInfo
    return data.computeValue(self)
  File "/home/felipe/.local/lib/python3.6/site-packages/ee/data.py", line 703, in computeValue
    prettyPrint=False))['result']
  File "/home/felipe/.local/lib/python3.6/site-packages/ee/data.py", line 340, in _execute_cloud_call
    raise _translate_cloud_exception(e)
ee.ee_exception.EEException: Image.sampleRectangle: Fully masked pixels / pixels outside of the image footprint when sampling band 'soilDiff' with no default value set. Note that calling sampleRectangle() on an image after ee.Image.clip() may result in a sampling bounding box outside the geometry passed to clip().

my exact code is the one below:

def imageToArray(img, geomet):
    band_arrs = ee.Image(img).sampleRectangle(ee.Geometry(geomet))
    print(band_arrs.getInfo())
    band_arr_b4 = band_arrs.get('eviDiff')
    band_arr_b5 = band_arrs.get('soilDiff')
    band_arr_b6 = band_arrs.get('mask')

    # Transfer the arrays from server to client and cast as np array.
    np_arr_b4 = np.array(band_arr_b4.getInfo())
    np_arr_b5 = np.array(band_arr_b5.getInfo())
    np_arr_b6 = np.array(band_arr_b6.getInfo())

    print(np_arr_b4.shape)
    print(np_arr_b5.shape)
    print(np_arr_b6.shape)

    np_arr_b4 = np.expand_dims(np_arr_b4, 2)
    np_arr_b5 = np.expand_dims(np_arr_b5, 2)
    np_arr_b6 = np.expand_dims(np_arr_b6, 2)

    rgb_img = np.concatenate((np_arr_b6, np_arr_b5, np_arr_b4), 2)
    print(rgb_img.shape)

    return rgb_img

Does anyone know how can I fix this?


texts: Creating buffers of a specific size and shape

texts:
I was recently asked to create a very odd "buffer" around several line features.  I'm not sure if this is innately possible with an ArcGIS Advanced license or not and therefore am willing to explore other solutions such as QGIS.  I am proficient in ArcPy, so please include any ArcGIS scripting solutions as well.

I have been given a set of lines that have a field called "Acreage".  The ultimate goal would be to generate polygons around each line that have the area of the "Acreage" field.  So, when considering the available data...  I have the length of each line and the area that each final polygon should be.  This alone doesn't seem like it should be to difficult to come up with a formula to calculate the buffer distance, especially if we were to assume each line is straight (though they are not straight, this might be the best I could do).  However, to make matters even more complicated those who asked for this would like to see the final result as "bent" rectangles.  The image below shows a mock-up of what they would like to see.

In the Image the red line is the starting data and the green "buffer" is what needs to be generated from each line.

The key elements are that the "buffer" needs to have the designated area(one of the line attributes) and also needs to have 90 degree corners.

If getting the corners is much more difficult or not possible.  I would also be happy with a solution that just creates a "normal" rounded buffer that contains the right area.


texts: Annotation boxes not drawing correctly in ArcMap 10.0 SP2

texts: I have about 400 annotation features in ArcMap 10.  The style is a simple text box with an arrow.  The annotation draws correctly when viewed up close on screen but the arrow detaches from the text box when printing or when zoomed out to view the entire 30x42" page.  I switched from labels to annotation to avoid problems like this.  Any ideas on fixing the printing issue?


texts: Seeking free high resolution imagery for Las Cruces, New Mexico area

texts:
I am looking for a free download Imagery that is a very high resolution and the most recently like the past 3 months.

I'm doing to produce a map where our new building being move into next year.

Is it possible for me to find one online?


texts: Carto DB: Scalable points over time using Torque in Carto DB

texts:
I am trying to visualize a data set of geo-points that have a change over time on an attribute. My goal is to have an animation that each point changes its scale over time.

I read the documentation but I couldn't find a solution.

The code that I tried was:

Map {
-torque-frame-count:10;
-torque-animation-duration:10;
-torque-time-attribute:"time";
-torque-aggregation-function:"count(were_here)";
-torque-resolution:1;
-torque-data-aggregation:linear;
}

#vivlio{
  comp-op: lighter;
  marker-fill-opacity: 0.9;
  marker-line-color: #FFF;
  marker-line-width: 1.5;
  marker-line-opacity: 1;
  marker-type: ellipse;
  marker-width: 6;
  marker-fill: #FF9900;
}
#vivlio[value>150] {
  marker-width:50;

}
#vivlio[value > 500] {
  marker-width:100;

}
#vivlio[value > 600] {
  marker-width:200;
}

An example of my data set is like this:

For a geopoint a

(Obj_id,Time,were_here)
(a,1,160)(a,2,600)(a,3,150)

In my animation I can't see scale difference over time.
Can anyone help me please?


texts:
I have an arcpy script which is bringing layers within a workspace into the table of contents, checking whether they contain features within the "current" extent of the mxd, and removing them if they are not. my script is as follows:

mxd = arcpy.mapping.MapDocument ("CURRENT")
arcpy.env.workspace = 'Database Connections\database.sde'
arcpy.env.overwriteOutput = True
df = arcpy.mapping.ListDataFrames(mxd, "*")[0]
arcpy.env.extent = df.extent

dsList = arcpy.ListDatasets("", "Feature")

for ds in dsList:
    fcList = arcpy.ListFeatureClasses("*","",ds)
    for FC in fcList:
        try:
            addLayer = arcpy.mapping.Layer(FC)
            addLayer.name = FC
            addLayer.visible = False
            arcpy.mapping.AddLayer(df,addLayer,"AUTO_ARRANGE")

            for lyr in arcpy.mapping.ListLayers(mxd,""):
                cnt = int(arcpy.GetCount_management(lyr).getoutput(0))
                if cnt<1:
                    arcpy.mapping.RemoveLayer(df,lyr)
                del addLayer
        except:
            print "problem with: " + str(FC)
            pass
            del FC
        del fcList

print 'finished'
arcpy.RefrechTOC()
arcpy.RefreshActiveView()
arcpy.AddMessage("Appropriate Feature Layers Added")
mxd.save
del mxd

This script works, but the trouble is, the database that I am working with has become larger than it used to be, so much so that this tool takes too long to run now for it to be useful.

Question: is there a way to call arcpy.mapping.RemoveLayer() outside of the 'for' loop of arcpy.mapping.ListLayers() in order to cut the size of the list the loop checks to 1?


texts: Leaflet activearea interferes with esri-lealet imagelayer

texts:
I am creating a leaflet application that uses various data sources / map layers.  I found a great source for global ground cover from esri.  It is an imageLayer.  Using an esri-leaflet imageLayer, I am able to add this layer to my map.  My application will likely need to draw data from various raster layers like this one.  I am also able to load an urban-wildland boundary layer as an esri-leaflet featureLayer.

I am also using leaflet-active-area, which is a plugin that shifts the center and bounds of the map to the div of your choosing.  I love this plugin for creating a better UI experience.

The featureLayer works just fine in conjunction with the activeArea.  However, the imageLayer is having problems.  It seems that the imageLayer is calculating what section of the layer should be currently in the map based on the original map bounds, but then squishing the image itself into the activeArea div.  Check it out:

Working demo of the problem

I made a codesandbox that shows this problem.  There is a 'toggle activearea' button there, so you can see how the imageLayer works properly when not using the activearea, but has this issue when using it.

What might be going wrong here?  What is it about an imageLayer (or its underlying L.ImageOverlay) that it knows to sample the correct bounds from the entire map, but then it squishes it into the activearea, instead of the total area?  I tried using other types of esri layers (L.esri.tiledMapLayer, L.esri.dynamicMapLayer, etc.) but these don't work as they are not the correct functions to recieve and apply this layer source.  I don't want to give up my activearea, but I'm going to need to apply and ready from raster layers in my application.   Ideas?

Thanks for reading.

Edit: altering codesandbox to use public esri imagemaplayer

I had been using an auth-requires esri imagemaplayer in my initial example, which had my personal client_id and client_secret exposed.  I changed it to use a public layer for esri so people can still view the problem.  Also,

Here is a codesandbox with the patch / fix js file

But thanks to TomazicM for getting to the crux of the matter.


texts: Use core code of QGIS in my own Python scripts

texts:
I would like to know if I can use the code of the default functions that come in QGIS, for example the Intersection function.
When you click on the help button it redirects you to the user guide where comes de description of the function and there is a section "Python code" where comes a small code for that specific function.
In almost all the vectors geoprocessing the library "processing" appears. But I could not install it from anywhere.
Does anyone know if it is possible use this code and library in my own Python script?


texts: How to perform GRASS v.rast.stats via QGIS

texts:
I have been using qgis for the last few months. I got interested in using GRASS tools through QGIS. I wanted to compute zonal statistics for a raster surface for a given polygon layer. I loaded the raster and the vector successfully into QGIS and GRASS. Then I used v.rast.stats Tool to compute zonal statistics and the program runs without any error. Here is a screenshot of the analysis:

But I can't see the results in my vector map. When I open the attribute table, no new column has been added. Can anyone advice me on how to get a new column with the computed value?

Thanks

Zaman


texts:
I've been trying to create a map and have it's initial zoom and center position be the query result's coordinates (a polygon, in this case) but it just refuses to zoom in, both using zoom in the constructor and centerAndZoom(). Centering, however, works just fine, it seems.

    _onBtnPrintClicked: function() {
                console.log("Begin query test");
                var queryUrl = "https://portalURL/arcgis/rest/services/FEATURE/FEATURE/FeatureServer/0";
                var queryTask = new QueryTask(queryUrl);
                var query = new Query();
                query.returnGeometry = true;
                query.outFields = ["*"];
                query.where = "OBJECTID = 53";
                console.log("Running execute");

                queryTask.execute(query).then(function(results) {
                    newExtent = new Extent(results.features[0].geometry.getExtent());
                    fichaMap = new Map("mapDiv", {
                        extent: newExtent,
                        showLabels: true,
                        logo: false,
                        zoom: 20,
                        slider: false
                    });
                }

Using Extent, like this, i get it to be centered and zoomed, but it is not the desired result as i would like it zoomed out a bit more (and not have the polygon in question fit the entire map's div element.
Now, i've already tried using center and zoom in the constructor and as well as using centerAndZoom and neither method works apart from centering.

What am I doing wrong?


texts: Identifying the greatest value point within each polygon

texts:
I have a series of polygons (4,000) and a series of points that exist on the perimeter of each polygon (3 to 5). The points refer to the elevation value which has previously been extracted from a DEM raster. The points exist at the minimum bounding locations of the polygons, whereby the Minimum Bounding Geometry Tool has been used to identify the outline of the polygons in a rectangular fashion and where each polygon intercepts its minimum bounding rectangle a point exists.
The polygons do not touch anywhere.
An example image of the points and polygons is shown below:

I wish to identify the point on the perimeter of each polygon that has the greatest elevation value. Any ideas how to do this in ArcMap 10.4?


texts: How do I compute v.kernel maps in less than 16 hours?

texts:
I am using the GRASS GIS 6.4.2 bundled with QGIS 1.8.0, using a PostGIS dataset exported via v.in.org.qgis from QGIS, to create kernel density maps.

The layer has 738,288 points.

My problem is the kernel density map takes about 16 hours to compute.

Here's parameters reported by v.kernel as it's executing:

STDDEV: 1000.000000
RES: 111.419043 ROWS: 458   COLS: 447

Writing output raster map using smooth parameter=1000.000000.

Normalising factor=6482635.018778.

I am running this on a Windows 7 x64 laptop with 8 GB RAM and an Intel Core i7 Q720 1.6 GHz with 4 physical cores. I notice that it's not multithreaded, only using 1 core.

Is there any way to speed this up without harming the quality of the output?


texts: Google maps and Ordnance Survey derived data licence restrictions

texts:
Current Status: You can't use data derived from restricted licence OS data with Google Maps.

There's been an on-going question about whether it is legal to display data derived from Ordnance Survey licenced mapping via google maps.

By derived data, I mean geometries which have been copied exactly from, or are composed of parts from, another dataset.

The blog post here: http://blog.ordnancesurvey.co.uk/2010/08/what-the-psma-really-means/ outlines the issue.

We've spoken to Google and they say that the OS have recently relaxed their licencing, but haven't been able to get any comment from the OS.  The last post on the OS blog is from November 2012.

Does anyone have any ideas as to how the land lies here?  Particularly if anyone has heard anything from the OS on this subject since November 2012.

UPDATE:
16th April 2013: We've got an account manager at the OS now and a collegue has just spoken to him.  Long story short: nothing has changed since November 2012 - any data derived from OS licenced mapping which is not specifically "free" (such as OS OpenData) may not be displayed in Google maps.  Our account manager didn't know when the OS was going to come back with an answer.  He said he'd chase it up.  *starts holding breath*

UPDATE:
2nd July 2013: Chasing our account manager, they reported the following (loosely quoted, do not base decisions solely on this) "The OS's position is unchanged and is not being pursued further. It is our (the data user's) responsibility to check with Google, or any other API provider, to ensure that OS rights are preserved. Rights to derived data cannot be given over to Google, as their T&C's require".  Now, talking to Google, they don't think their T&C's violate the OS licence.  I'm not sure this is too useful though.  Our account manager has promised to send written clarification, which I will post here: a) If it ever arrives and b) if it isn't confidential.  The OS blog hasn't been updated yet.

UPDATE:
2nd July 2013: Just had an email from our account manager which contradicts what we were told earlier: "...have had a chat with our lawyers who confirmed that we are in discussion with Google around the compatibility of their terms with ours.  Currently, the position is that Ordnance Survey map data cannot be overlaid on top of Google Maps though we will not have an “official document/wording” available until the current discussions are concluded.".  So there!

UPDATE:
6th September 2013: An email to a colleague from the boss of our account manager includes the line: "...Currently our data cannot be overlayed on top of Google Maps due to their rules around IP ownership".  I don't know if this is now their official legal stance or this guy's opinion.


texts: Spatial Precision for Editing in Multiple GIS Clients

texts:
I have something I'd like to discuss, I'll try to explain it first, why it's important, and I'll explain how ArcGIS handles the issue, which I think is a good way to handle it. Finally I'd like to know what you think.

The issue is this:

Many GIS editing software packages can edit many data formats. Each data format handles precision differently, for example, SHP format stores a double precision coordinate. SQL Server does the same, a 64 bit double. PostGIS seems to do the same. MapInfo's .TAB format apparently uses a 21 bit int to store a point as a linear interpolation between bounding coordinates. There will obviously be others. Editing software may have internal data structures whilst carrying out operations, and will then write to the end data store e.g., MapInfo linked tables store a value on a (for us, a 1 mm) grid, and push to the database, which stores it as a double precision float.

If you edit in mixed environments, then you can start getting inconsistencies where a feature with coordinates that were created off a grid (with doubel precision) are then edited in software that uses a grid (in mm precision). For example, I create a layer in SQL Server, add features in QGIS, then edit one of the vertices of one feature in MapInfo, that feature now gets all its coordinates snapped to the 1 mm grid, which means it doesn't align to neighboring features any more (the errors in my test were ~0.2 mm).

ArcGIS handles precision when using a database layer by keeping track of xy_resolution (see here for more on xy_resolution and xy_tolerance) explicitly. This seems the best way of handling it, since the coordinate precision is explicitly controlled (the default is 1 mm). When editing in MapInfo with carefully specified bounds you snap to a 1 mm grid implicitly, so this is sound too.

My question is:

How do you or your organisation handle editing in order to maintain topology, keeping in mind different sources with different precisions?

And as a secondary question, should editing tools, such as QGIS, include an option to use a defined coordinate precision or should we do it on data format, with a database trigger or MapInfo's bounds?

Regards,

Alex

--Further Details follow:

I have two images, one of the shape I'm editing, one of the error introduced if editing with different precisions. I also have numbers showing what is being written to SQL Server when you move the central point in two different softwares.

QGIS

POLYGON ((520191.61484399455 5259827.5143016046, 521146.0563152277 5259917.96005114, 521015.6907808886 5260649.8136086222, 520329.84600000043 5260729.1830000039, 520191.61484399455 5259827.5143016046))
POLYGON ((521146.0563152277 5259917.96005114, 522006.03708566393 5259775.4479602445, 522021.9278674855 5260443.1099995337, 521015.6907808886 5260649.8136086222, 521146.0563152277 5259917.96005114))
POLYGON ((521938.09510281781 5258959.9356564963, 522006.03708566393 5259775.4479602445, 521146.0563152277 5259917.96005114, 521194.71867216442 5258970.2782813907, 521938.09510281781 5258959.9356564963))
POLYGON ((521194.71867216442 5258970.2782813907, 521146.0563152277 5259917.96005114, 520191.61484399455 5259827.5143016046, 520191.61484399455 5259827.5143016046, 520287.61066861876 5259041.6976504158, 521194.71867216442 5258970.2782813907))

MapInfo

POLYGON ((520191.615 5259827.514, 520842.732 5259437.788, 521015.691 5260649.814, 520329.846 5260729.183, 520191.615 5259827.514))
POLYGON ((520842.732 5259437.788, 522006.037 5259775.448, 522021.928 5260443.11, 521015.691 5260649.814, 520842.732 5259437.788))
POLYGON ((521938.095 5258959.936, 522006.037 5259775.448, 520842.732 5259437.788, 521194.719 5258970.278, 521938.095 5258959.936))
POLYGON ((521194.719 5258970.278, 520842.732 5259437.788, 520191.615 5259827.514, 520287.611 5259041.698, 521194.719 5258970.278))

Polygon that was editied, the central point was moved.

The resulting difference after editing in MapInfo, the dataset is no longer topologically correct.


texts: How to incorporate drop-down list within cells of shapefile attribute table using ArcGIS Desktop?

texts:
A question has been brought to my attention. When working with a shapefile within either 10 or 10.1. Is it possible to create a python script or some sort of programming language to accomplish this task;

An attribute table consists of 3 columns with multiple rows. Is there a way to set a list for each column that may be used when clicking within each cell. I know this can be done within an edit session with domain and subtypes but can it be done just within the Attribute table itself. By way of scripting.

Reason, this would provide a streamline process for individuals who do not know much about ArcGIS and would save the process of an edit session and potential of misspellings and what not. The easy fix for a GIS person would be to just use Domains and Subtypes. But the question is, can it be done differently?


texts: NetworkX: Saving Shortest path of OpenStreetMap data as an OGR friendly format

texts:
I am using this script to import OSM data into a GIS based model.

This has to run independent from any software package such as QGIS, PostGIS, or ArcGIS.

I then create a shortest path between two nodes as follow:

shortest_path = nx.shortest_path(G, source=u'879835512', target=u'879912131')
result_graph = G.subgraph(shortest_path)

I would like to visualize (kml?) the results of the shortest path and store them (shp?).

I have a feeling that I need to use the write_shp but I think that only works if you have already used read_shp which I haven't since I'm importing directly from OSM.

EDIT 1:

Based on @gene 's comment, I added the following code:

for e in result_graph.edges():
    p0x = result_graph.node[e[0]]['lat']
    p0y = result_graph.node[e[0]]['lon']
    p1x = result_graph.node[e[1]]['lat']
    p1y = result_graph.node[e[1]]['lon']
    edges.append(LineString([(p0x, p0y), (p1x, p1y)]))

nodes = []
for n in result_graph.nodes():
    px = result_graph.node[e[0]]['lat']
    py = result_graph.node[e[0]]['lon']
    nodes.append((px,py))

G2 = nx.DiGraph()

G2.add_nodes_from(nodes)

I don't think I'm on the right track yet. as it seems I need to maintain the coordinates somehow with my nodes/links


texts: How can I calculate times with gps data using random points on the track?

texts:
I have GPS data from a car that is circulating all the time in a closed circuit. I need to put a few points along the route, about 3 or 4 and calculate the time it takes to get between points. These points would be 3 or 4 points, no more. And also for the whole lap time.

How can I calculate the time between those points at random positions and the calculated time to make a full lap?  Maybe I need to use Excel for example?

Someone recommended me to create a Network with circuit lines that follow as a road network. But I really do not understand the concept.


texts:
How can I validate the geometry in Openlayers?

I need to check  "MULTILINESTRING()" (which is empty) is valid/not valid. In postGIS, for Geometry validation we are using "isValid()" function. Like this, is there any options in Openlayers?


texts: Segmentation in OTB gives blank rectangle as output

texts: I am trying to segment a raster image to get a vector output. I am using OTB 6.6.1 in QGIS 3.4. In the segmentation window, I selected the input raster which is a 3 band (true colour) sentinel 2 image and chose the meanshift under segmentation algorithm and chose an output file name with .shp as output. I left all the other options to default. I am getting the output as a blank rectangle. I am getting a perfectly segmented image if I choose watershed instead of meanshift in segmentation algorithm. Any suggestions ?


texts: Empty layers while processing 'extract by location' in Python

texts:
I have a folder with subdirectories containing shapefiles. I iterate in those folders to create a dictionary named liste_couches (['name of layer': complete path to shapefile']).

The code works fine when I directly call it, it returns a layer called 'fichiertraite' with the good extraction:

layer1 = QgsVectorLayer(liste_couches['geobase'], 'SamGamegie', 'ogr') #layer to be clip # liste_couches is a dictionary ['name of layer': complete path to shapefile']

tmp = processing.runalg('qgis:extractbylocation', layer1, lyTampon, [u'intersects'],0,None) #lyTampon is my intersect layer
lyc = QgsVectorLayer(tmp["OUTPUT"].decode('utf8').encode(enco), 'fichiertraite', 'ogr')
QgsMapLayerRegistry.instance().addMapLayer(lyc)

But when I nest it like this in a loop, it returns empty layers:

for couche, chemin in liste_couches.iteritems():

    lyCouche = QgsVectorLayer(chemin, couche, 'ogr')

    QgsMapLayerRegistry.instance().addMapLayer(lyCouche)
    tmp = processing.runalg("qgis:extractbylocation", lyCouche, lyTampon, [u'intersect'],0,None)
    lyc = QgsVectorLayer(tmp["OUTPUT"].decode('utf8').encode(enco), couche + '_tr', 'ogr')
    QgsMapLayerRegistry.instance().addMapLayer(lyc)


texts: rgdal loads PROJ.4 version 4.9.1 when the updated 4.9.2 is installed

texts:
I have been having problems with vectorial spatial data for some time, mostly important to plot SpatialPolygonDataFrames (I asked here: Slow plotting SpatialPolygonDataFrames on some systems). This happens to me in a Macos laptop (MacBook Pro (Retina, 13-inch, Late 2013) with macOS Sierra 10.12.4 (16E195).

Trying to solve it, I updated all libraries related to spatial data, among them GDAL and PROJ. After reinstalling rdgal, I have noticed that when I load the library, the wrong version of PROJ.4 is loaded:

> library(rgdal)
rgdal: version: 1.2-7, (SVN revision 660)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.2, released 2016/10/24
 Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
 Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
 Linking to sp version: 1.2-4

I am sure I have installed the version 4.9.2, and it must be properly set in PATH because if I run the command proj on the terminal I get:

MacBook-Pro:~ javierfajnolla$ proj
Rel. 4.9.2, 08 September 2015
usage: proj [ -beEfiIlormsStTvVwW [args] ] [ +opts[=arg] ] [ files ]

The version 4.9.1 is known to have some problems (although I totally do not know if my problems are at all related to those... I'm just desperate for a solution), so I think I should update it. But I have failed trying to make rgdal read the right version of it.

Does anyone know a way to do this?

EDIT: I still do not have a way to solve this issue, but I have some additional information:

If the plots are done in R, not RStudio, it is faster. Not as fast as in other OS, but still better. It seems that the plotting device built in RStudio is quite slow.
There is a workaround related with the previous point. It is possible to use a different plotting device that RStudio's, and that makes it much faster. You can open pdf() or png(), make the plot and close the connection dev.off() and you will get a file with the plot inside. Not ideal, but certainly working.


texts: Why is this field calculator expression evaluated incorrectly when using Python syntax?

texts:
I have a polygon featureclass containing an integer field 'gridcode'. I've added 2 new float fields VB_calc and Python_calc, both of which should be equal to gridcode / 1,000,000.

If I calculate the field using the Field Calculator with the VB Script [fieldname] syntax, the result is as expected.

If I use the Python !fieldname! syntax, the results are all zero.

What am I doing wrong here? I entered those expressions by double-clicking the gridcode fieldname in the calculator. (My goal is to use this within a stand-alone Python script but I'm debugging it here since the script isn't working either.) Both fields are definitely of type FLOAT.


texts:
Problem:

I'm trying to loop through all my layers looking for layer named "CADAnnotation".
If the layer exists then remove the layer from the mxd

Notes:

Running from Stand-Alone script (ie. NOT within Arcmap)
"CADAnnotation" Data Type is a CAD Annotation Feature Class
"CADAnnotation" is NOT in a geodatabase, it's created from a AutoCAD .dwg
"CADAnnotation" is within a Group Layer named "ACAD"
If the group layer "ACAD" can be deleted that also removes "CADAnnotation" that would be great.

Code thus far:

for item in mxds:
    print (item)
    mxd = arcpy.mapping.MapDocument(item)
    df=arcpy.mapping.ListDataFrames(mxd,"Project Area")[0]
    for lyr in arcpy.mapping.ListLayers(mxd, "*",df):
        if lyr.name == "CADAnnotation":
            print(lyr.dataSource)
            arcpy.Delete_management("CADAnnotation")
            print("Layer Deleted")
        else:
            pass

Notes on Code:

I can find the layer no problem
the line arcpy.Delete_management("CADAnnotation") does not work throws an error.

Question:

How to I remove "CADAnnotation" and/or "ACAD" group layer?


texts:
I know how to set projection and georeference bounds of scanned TIFF file, by using gdal, but I don't know how to do the same for vector file.

For example I georeferenced TIFF file by using this command:

gdal_translate -of GTiff -a_ullr x1 y1 x2 y2 -a_srs "7041.prj" input.tif output.tif

where 7041.prj is Gauss Kruger projection file.

I used CorelDraw to do the tracing, and both raster and DXF output share same bounds - vector file has boundbox same as raster image used for tracing. I wanted they overlay automatically when loaded in GIS.

For vector file, with ogr2ogr I can convert DXF to SHP and set projection, but I can't figure how to set georeference bounds as ogr2ogr doesn't accept -a_ullr argument, probably for some reason.

So my question is can I somehow set bounds for resulting SHP file?


texts:
I recently went through a lengthy process of installing GDAL and think that something may have gone awry as my IDLE is not opening. I cannot edit files in IDLE either. The only thing that works is running the IDLE.bat which then opens the Python Shell. I do have Python 2.7 installed via Arc (thus my installation for Python is in C:\Python27\ArcGIS10.4).

I've tried copying the tcl8.5 folder into the Lib folder. I've cleared the files in my .idlerc folder. Short of reinstalling Arc entirely, are there any workarounds? I can't seem to find much else that helps.As you can see in the photos, my .py files have also lost their default programs even though pythonw is set as the default program.

Can't seem to add a proper photos link, but the photos are at http://imgur.com/a/eO1Tr


texts:
I have a styling function like this:

var customStyleFunction = function( feature, resolution ) {
       var resultStyles = [];

       (other styles, line, point, etc)

       var featureText = new ol.style.Style({
           text: new ol.style.Text({
               text: feature.getProperties().label,
               scale : 1,
           })
       });

       resultStyles.push( featureText );
       return resultStyles;
}

When I have few polygons, All goes ok.

but when I try to add a vector layer with some more polygons, my map crashes:

the viewport will reduce every time I try to move or zoom the map. What I have done wrong?


texts:
I need a help with my homework.
I have created digital elevation model from the shapefile. In QGis(Grass plugin)

This is how far i can go. I have to find the longest slope oriented to North.
I know I have to use r.slope and r.aspect modules, but how can I use it?

Thanks!


texts: Calculating area of polygon considering relief in ArcGIS Desktop?

texts:
I work in a fire department and i'm mapping all the fires that occur in my area, creating polygons that represent burnt areas.
my problem is that the values are not accurate... for instance, arcgis calculates one of the polygons as 51000 square meters when, in reality, I know that the real area is 60000 square meters. the issue is that the relief is not being considered for the measurement and I need to have the most accurate possible values.

is there any way to calculate an area considering the relief in arcgis 9.3?


texts: Adding layer to open ArcMap document using ArcPy?

texts:
I have created a script. I'm run the script using ArcMap and a layer should be added to the open ArcMap document.

I have earlier used MXD.save() witch work fine when I'm running the script outside ArcMap (through an IDE).

How can I change the code below such that the layer is added to the open ArcMap document?

Code:

MXD = arcpy.mapping.MapDocument(r"C:\Users\mgv\Desktop\Koncept\Koncept.mxd")
DF = arcpy.mapping.ListDataFrames(MXD, "Layers")[0]
arcpy.MakeFeatureLayer_management(in_features=line + "_autodimensioned", out_layer=line + "_autodimensioned_layer")
layer = arcpy.mapping.Layer(line + "_autodimensioned_layer")
arcpy.mapping.AddLayer(DF, layer, add_position="AUTO_ARRANGE")

# Add symbology for the layer showing the pipes, symbology for the pipe_type_sup field
updateLayer = arcpy.mapping.ListLayers(MXD, "C:\Users\mgv\Desktop\Koncept\Network.mdb\Network\line_autodim_autodimensioned_layer", DF)[0]
sourceLayer = arcpy.mapping.Layer(r"C:/Users/mgv/Desktop/Koncept/Symbology_configuration_layer.lyr") #task: aendre
arcpy.mapping.UpdateLayer(DF, updateLayer, sourceLayer, True)

# Add label with the pipe type for the supply pipe (DN)
lcs = updateLayer.labelClasses
lc = lcs[0]
lc.expression = "[pipe_type_sup]"
updateLayer.showLabels = True

MXD.save() #Works if running in PyCharm / running the code outside ArcMap / ArcGIS for Desktop


texts: Writing GeoJSON stream into Shapefile using GeoTools?

texts:
I've been trying a few solutions to efficiently convert GeoJSON into a
shapefile without having to store all features in memory. I'm using
GeoTools 9.2.

The problem is not so much in how to stream the JSON but how to
efficiently write the features into the shapefile. I use
FeatureJSON#streamFeatureCollection to obtain an iterator. After some
googling, I found 3 different ways of writing a shapefile, namely:

Option 1: Repeatedly calling FeatureStore#addFeatures with a collection containing say 1000 features, within a transaction.

  ListFeatureCollection coll = new ListFeatureCollection(type, features);
  Transaction transaction = new DefaultTransaction("create");
  featureStore.setTransaction(transaction);
  try {
    featureStore.addFeatures(coll);
    transaction.commit();
  } catch (IOException e) {
    transaction.rollback();
    throw new IllegalStateException(
        "Could not write some features to shapefile. Aborting process", e);
  } finally {
    transaction.close();
  }

This option is extremely slow. By profiling a few runs, I noticed that
about 50% of CPU time is spent on the method
ContentFeatureStore#getWriterAppend, presumably in order to reach the
end of the file before each transaction commit.

Option 2: Obtaining an append writer directly from ShapefileDataStore, and
write 1000 features at a time within a transaction.

This option suffers from the same problems as number one.

Option 3: Obtaining a feature writer from ShapefileDataStore, and write one feature at a time using Transaction.AUTO_COMMIT.

 FeatureWriter<SimpleFeatureType, SimpleFeature> writer = shpDataStore
    .getFeatureWriter(shpDataStore.getTypeNames()[0],
        Transaction.AUTO_COMMIT);

 while (jsonIt.hasNext()) {

  SimpleFeature feature = jsonIt.next();
  SimpleFeature toWrite = writer.next();
  for (int i = 0; i < toWrite.getType().getAttributeCount(); i++) {
    String name = toWrite.getType().getDescriptor(i).getLocalName();
    toWrite.setAttribute(name, feature.getAttribute(name));
  }
  writer.write();
}
writer.close();

Option 3 is the fastest, but I feel there would a way of efficiently
adding a greater number of features at a time to the shapefile within
a transaction.

On the other hand, a previous comment in the GeoTools mailing list
noted:

  The above would work for mid-sized data transfers, for massive ones
  against databases it's better to adopt some sort of batching to avoid
  having a single transaction with one million inserts, e.g., insert
  1000, commit the transaction, insert another 1000, and so on. This
  would work better against databases and against WFS servers, but not
  against shapefiles, which instead work better with the massive
  insert... to each his own.

Does this mean that the most efficient way of writing to a shapefile
is having all features in memory, rather than being able to append
features?


texts: Deleting all lines not perpendicular to another line layer using QGIS?

texts:
I am using the latest QGIS 3.4.4

I have a layer 1 which has a bunch of lines in all directions. I have another layer  2 intersecting those lines.

I want to delete all the lines in layer 1 that do not cross perpendicularly to layer 2.

Data: You can download the two line layers here

This seems simple but I cant seem to get it right.


texts: Incorrect cell count with viewshed analysis

texts:
I am performing several viewshed analyses. As I wanted to observe how the viewshed changes when setting a different value for OffsetA I noticed something strange was going on with the cell count.

Firstly, I performed the analysis without an offset value. These were the results:

The attribute table states that the cell count is 208141.
Then I performed another analysis with an OffsetA of 2m, and I got these results:

Obviously, the viewshed is larger (which I expected it to be), but the cell count as mentioned in the attribute table is lower (113111).

Since I set invisible cells to be returned as NoData, I suspected it might had something to do with that. Therefore I performed the analysis again, but without selecting the option for invisible cells to be displayed as NoData.

This did not change anything and the results were exactly the same.

Am I doing something wrong or am I overlooking something?
I am using ArcMap 10.3.


texts: Is there GIS data for Electrical Cables under ground in US?

texts:
Where is GIS data for electrical cables under the ground in US ?
AABSyS seems to provide data, but I couldn't find out it after all.
The data might not be in public. In that case,
substitute to it, such as trans (under/above ground),
railway and so on, works for me.

Since I'm using QGIS, any data which can be imported into QGIS is preferable.


texts: What does this warning mean for shapely python

texts:
The program keep outputting:

/Users/opt/anaconda3/lib/python3.8/site-packages/geopandas/_compat.py:106: UserWarning: The Shapely GEOS version (3.8.0-CAPI-1.13.1 ) is incompatible with the GEOS version PyGEOS was compiled with (3.9.1-CAPI-1.14.2). Conversions between both will be slow.
warnings.warn(

when the program below is written:
   import numpy as np
   import pandas as pd
   import matplotlib.pyplot as plt
   import geopandas as gpd
   from descartes import PolygonPatch
   from shapely.geometry import Point Polygon

I have tried to solve this researching for hours but don't know what it means and when I try to do any program related to shapely the incorrect solution is outputted.


texts: Create a line from points using cartodb.js

texts:
I would like to draw a line between multiple points without creating a new dataset. Using this tutorial I've been able to draw lines between points in the CartoDB Editor:
http://docs.cartodb.com/tutorials/gps_track/

However, I would like to do the same using Cartodb.js. But I can not get this to work. When I have drawn the lines in the layer in CartoDB editor it's also working on cartodb.js, probably because the connected JSON layer URL then understands that it should work with lines. However when I clear all SQL in CartoDB editor and see the points again, I can not 'convert' it to lines using Cartodb.js.

This is the code I'm using. Is there a fix so I can create lines directly from the dataset which is containing points?


texts:
I'm just beginning to play with BoundlessGeo's GeoExplorer interface. I'm trying to figure out if I can customize it for use in-house. We basically need an intra-net type web mapping app to display various data layers (client project info, location, etc.). We need to do some custom queries, and there;s interest in a tabbed identify display.

Anyone know if it can be customized in such ways?

I'm finding the documentation kind of lacking, and this seems like a logical place to turn for answers.


texts:
I am updating construction works using QGIS. Every 6 weeks I update the map with new changes, but need to save each survey cycle as a new file. I have a map set up already with all the layers etc.

My issue is that I'm about to start a new survey cycle and need to save the changes I find. I have saved the project for August into a new file, but how do I export all the existing layers etc into a new folder so that when I make changes, it will be for new layer files instead of updating the layer files for the previous batches project?


texts:
I need to subset a global fishing dataset to just have the fishing in the Pacific. Right now I have a shapefile of the Pacific and I am using the sp method of subsetting (see below)

oceans <- readOGR(dsn = ".", layer = "PacificOcean")
gfw = readRDS(filename2)
gfw <- SpatialPointsDataFrame(coords = gfwcoords, data = gfwdata,
                              proj4string = CRS(proj4string(oceans)))
gfwmask<- gfw[oceans, ]

This works just fine, BUT, I need to do it hundreds of times and it takes several minutes each time. Does anyone have suggestions for a faster method to subset a dataframe by a shapefile?


texts: ArcGIS is not recognizing my network connections. What do I do?

texts:
So, I've noticed this weird thing when I run ArcMap as administrator.  It shows my folder connections as not connected. The program does not seem to recognize that I have any network drives at all.

Has anyone else ever had this problem, and how do I fix it?

Running 10.1 on Windows 8.1 64 bit

Update:  When I go to make a new folder connection when running as administrator, if I try to click on network locations the program errors out and crashes.


texts: Calculating polygon attribute based on intersection with grid in QGIS

texts:
I'm trying to use QGIS to profiling polygons in a layer. There is also another grid layer, each grid contains its own demographics information.

Assume the demographics are evenly distributed within the same grid. Therefore, I'm deriving the value for any demographic attribute of an outlet by summing the sections in different grids weighted by the area.
and the calculation for area A, for instance, population, would be:

I know I can use the clip to have the intersects layer, use the field calculator to get the areas and then calculate the population for A1 and A3 separately.

How could I aggregate the sum population and add the population back to area A in the original polygon layer?


texts: Getting the Intersection of lines not at same location

texts:
I have got area with line strings that are not quite similar to each other means that at some places they are not exactly at same position. This can be seen in the picture. I am making a buffer on linestring Area1 and then comparing it with the other linestring Area2. I am using ST_Intersection to get the intersection of buffer and linestring. As both areas are not located on each other so I have to increase the size of the buffer. There are also some extra linestrings in Area2 which I dont want. I only want the geometry of Area2 not Area1. My problem is that once I increase the size of buffer to get the intersection of both areas the extra lines in Area2 also comes in the geometry. I either want to remove these extra linestring from geometry or they should be compared in such a way that extra lines do not come from area2.

The first picture shows how area2 and area1 differs.

The next pictures shows that If I increase the size of buffer it will include extraline which I dont want. The orange line is from Area2 and purple buffer is Area1

The query I am using is this in which I am getting the extra parts from AREA2 which I dont want.

Create Table Mytable As select
ST_MULTI(ST_Intersection(Area2,ST_Buffer(SetSrid(Area1.geometry,3857),8)))  AS Geometry
from Area2 ,Area1  where
ST_Intersects(Area2.geometry,ST_Buffer(SetSrid(Area1.geometry,3857),6))

I have tried everything but could not find the way.

Can Someone tell its not working? I have already tried ST_SNAP and ST_SNAPTOGRID but it is not working


texts:
I have a dataset with 2157335 entries, and I am performing IDW/Kriging to interpolate.

But the code is running very slow. I am not sure this is because of a large dataset or there is something wrong with my code.

grd <- Sobj_SpatialGrid(ndvi_pts, maxDim=100)$SG
response <- unclass(ndvi_pts@data)[[1]]
ndvi_idw <- idw(response~1, ndvi_pts, newdata=grd)

or for kriging,

kriging_result=autoKrige(response~1, ndvi_pts, new_data=grd, nmax=7)

I cannot see anything wrong with it and R doesn't show any error. It just keeps running with no sign of stopping.


texts:
I'll try to keep this simple.

Is there a way to make a polyline that runs along the edge of these polygons. It is cadastrals that assembly the coastline and I need to generate a line that works as the outer boundary.


texts: ArcGIS Rest API returns No Results When Should Intersect

texts:
I can download statewide NY senate districts using ...

https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/tigerWMS_Current/MapServer/56/query?where=STATE+%3D+36&text=&objectIds=&time=&geometry=&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=*&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&resultOffset=&resultRecordCount=&queryByDistance=&returnExtentsOnly=false&datumTransformation=&parameterValues=&rangeValues=&f=geojson

When I query lat/lng against the api... I get zero results.

https://tigerweb.geo.census.gov/arcgis/rest/services/TIGERweb/tigerWMS_Current/MapServer/56/query?where=STATE+%3D+36&text=&objectIds=&time=&geometry=-73.831676%2C42.768687&geometryType=esriGeometryPoint&inSR=&spatialRel=esriSpatialRelIntersects&relationParam=&outFields=*&returnGeometry=true&returnTrueCurves=false&maxAllowableOffset=&geometryPrecision=&outSR=&returnIdsOnly=false&returnCountOnly=false&orderByFields=&groupByFieldsForStatistics=&outStatistics=&returnZ=false&returnM=false&gdbVersion=&returnDistinctValues=false&resultOffset=&resultRecordCount=&queryByDistance=&returnExtentsOnly=false&datumTransformation=&parameterValues=&rangeValues=&f=geojson

But when I map the statewide dataset and map the lat/lng I can see I should have an intersect.  Can someone tell me why this intersect returns no results when the point falls within zone on the map???


texts:
While I'm working in post gis sql command prompt its mode changing automatically.

I mean that the general form of command prompt as follows

template_postgis=#(we type query's here)

While working some time, when a query given its changed as

template_postgis<#

if any key pressed its unchanged again comes as

template_postgis<#

how can I change this to normal form like template_postgis=#?

psql (8.4.5)
WARNING: Console code page (437) differs from Windows code page (1252)
         8-bit characters might not work correctly. See psql reference
         page "Notes for Windows users" for details.
Type "help" for help.

template_postgis_20=# create table cities (id int, name text, state_name text);
CREATE TABLE
template_postgis_20=# select addgeometrycolumn ('', 'cities')
template_postgis_20-# select addgeometrycolumn ('', 'cities', 'latlon', 4326, 'M
ULTIPOINT', 2);
ERROR:  syntax error at or near "select"
LINE 2: select addgeometrycolumn ('', 'cities', 'latlon', 4326, 'MUL...
        ^
template_postgis_20=# SELECT addgeometrycolumn ('', 'cities', 'latlon', 4326, 'M
ULTIPOINT', 2);

                   addgeometrycolumn

public.cities.latlon SRID:4326 TYPE:MULTIPOINT DIMS:2
(1 row)

template_postgis_20(# INSERT INTO cities (id, name, state_name, latlon) VALUES (
'1', 'city', 'andhra', GeometryFromText('MULTIPOINT(78.4667 17.3667,82.2200 16.9
300,81.7833 16.9833,79.3515 13.6797,83.0000 17.0000,80.6200 16.5200)', 4326);
template_postgis_20(# INSERT INTO cities (id, name, state_name, latlon) VALUES (
'1', 'city', 'andhra', GeometryFromText('MULTIPOINT(78.4667 17.3667,82.2200 16.9
300,81.7833 16.9833,79.3515 13.6797,83.0000 17.0000,80.6200 16.5200)', 4326);
template_postgis_20(# INSERT INTO cities (id, name, state_name, latlon) VALUES (
'1', 'city', 'andhra', GeometryFromText('MULTIPOINT(78.4667 17.3667,82.2200 16.9
300,81.7833 16.9833,79.3515 13.6797,83.0000 17.0000,80.6200 16.5200)', 4326);
template_postgis_20(# INSERT INTO cities (id, name, state_name, latlon) VALUES (
'1', 'city', 'andhra', GeometryFromText('MULTIPOINT(78.4667 17.3667,82.2200 16.9
300,81.7833 16.9833,79.3515 13.6797,83.0000 17.0000,80.6200 16.5200)', 4326);

template_postgis_20(#
template_postgis_20(# INSERT INTO cities (id, name, state_name, latlon) VALUES (
'1', 'city', 'andhra', GeometryFromText('MULTIPOINT(78.4667 17.3667,82.2200 16.9
300,81.7833 16.9833,79.3515 13.6797,83.0000 17.0000,80.6200 16.5200)', 4326);
template_postgis_20(# fgfgfd
template_postgis_20(# fgfdgfd
template_postgis_20(# gfdg
template_postgis_20(# re
template_postgis_20(# tre
template_postgis_20(# tr
template_postgis_20(# gfgretrytojyhftipg ccnb ds
template_postgis_20(# fgfdgfd
template_postgis_20(# fgfdgfd
template_postgis_20(# 222
template_postgis_20(# @@@@@@@@@@@@@@@@
template_postgis_20(#
template_postgis_20(#
template_postgis_20(#


texts: Plot GPX on topographic map with no labels

texts: I am trying to plot a path from a GPX file on a topographic map with as few labels as possible. I've played around with a handful of both online and desktop solutions, but it seems that I either don't have enough control over what labels do and don't show up on a given map, or the labels are native to the tiles themselves. The final product of this will be a map that I'll print on letter-sized paper so that we can show it to people as part of a research project. We're willing to spend some money for a piece of software if it will solve this for us in a relatively easy manner. Is there a tool, online or otherwise, that will allow us to do this?


texts: How to show GPS latitude/longitude in input in Android application (webview)

texts:
I develop an Android application with jQuery Mobile 1.9, Cordova 2.4 ... I want to display the coordinate of GPS (lon,lat) in input like in the pictures, I'm tired with this code but it can't display the coordinates. If you know where error is help me please.

.


texts:
I just installed QGIS 3.0 EA.

I wonder where the Spatial query is?

There is a ”count points in polygon” in the toolbox but it can’t count point within a polygon.


texts: Addition of line of latitude to a coordinate system

texts: I need to know how to add a line of latitude to a graticule. I want to add a line of 66 deg. 34 minutes to north pole orthographic graticule. I am using ArcGIS


texts: How to find data about all individual residential properties in a neighborhood and map them out parcel by parcel?

texts: I am trying to make a map about residential properties in 3 neighborhoods in Oakland, California. I already have the shapefile for all the parcels and blocks, but I need to find data that indicate which and which specifically is a residential property, and highlight them all to see the distribution of residential propoerties.


texts:
I tried to produce a smooth line animation with below code. How can I do this?

I used OpenLayers with Angular. It works as a drawing line but I want more smooth movement something like this (http://openlayers.org/en/latest/examples/feature-move-animation.html). I didn't modify my code because this example have JS but I am working on TS. Besides I used random data which is coming from Python client by websockets to node.js server.

I just want to flowing marker, but in my code marker has stopped and then appear again. I need consistent marker and it should not disappear.

map.component.ts

import { Component, OnInit, Input, Output, SimpleChanges, EventEmitter } from '@angular/core';
import * as ol from 'openlayers';
import { AppStateService } from '../state/appstate.service';
import { DataService } from 'src/core/data.service';
import {Subscription, timer} from 'rxjs';
const INITIAL_OPACITY = 1;
const DIMMED_OPACITY = 0.3;
var longval = 0;
var latval = 0;
@Component({
  selector: 'map',
  template: `<div id="map"></div>`,
  styleUrls: ['map.component.css']
})
export class MapComponent implements OnInit {
  private map: ol.Map;
  stockQuote: number;
  sub: Subscription;
  constructor(private appStateService: AppStateService, private dataService: DataService) {
  }
  ngOnInit() {
    this.map = new ol.Map({
      layers: [
        new ol.layer.Tile({ source: new ol.source.OSM(),opacity: INITIAL_OPACITY }),
      ],
      target: document.getElementById('map'),
      view: new ol.View({
        projection: 'EPSG:4326', // map type
        center: [32.7399307, 39.9066235],
        zoom: 6,
      })
    });
    this.appStateService.getDim().subscribe(
      x => {
        if (x) {
          this.map.getLayers().getArray()[0].setOpacity(INITIAL_OPACITY);
        } else {
          this.map.getLayers().getArray()[0].setOpacity(DIMMED_OPACITY);
        }
      }
    )
    this.sub = this.dataService.getQuotes()
      .subscribe(value => {
        var long = value['long'];
        var lat = value['lat'];
        if (longval === 0)
          longval = long;
        if (latval === 0)
          latval = lat;
        var source = new ol.source.Vector({});
        var layer = new ol.layer.Vector({ source: source });
        this.map.addLayer(layer);
        var marker = new ol.Feature({
          zoom: 8,
          geometry: new ol.geom.LineString([[lat, long], [latval,longval]])
        });
         var marker2 = new ol.Feature({
          zoom: 8,
          //geometry: new ol.geom.Point([lot, lang]) // dont worry about coordinate type 0,0 will be in west coast of africa
          geometry: new ol.geom.Point([lat, long])
          //geometry: new ol.geom.LineString([[lat, long], [latval,longval]])
        });
        source.addFeature(marker)
        function delay(ms:number){
          return new Promise( resolve => setTimeout(resolve, ms) );
        }
        (async () => {
            // Do something before delay
            source.addFeature(marker2)
            await delay(500);
            // Do something after
          source.removeFeature(marker2)
        })();
        longval = long;
        latval = lat;
      });
  }


texts: Exit searchCursor loop after first result

texts:
I am attempting to join 2 tables together, and the tools in arcmap, such as add field can't do what I want so I tried using python and arcpy. My approach is to use 2 nested loops and find the values in the 2nd table and update the first one. This is the code that I am using:

with arcpy.da.UpdateCursor(table1, ufield) as ucursor:
for urow in ucursor:
    print(urow[0])
    with arcpy.da.SearchCursor(table2, sfield, """"VALUE" = """ + str(urow[0])) as scursor:
        for srow in scursor:
            print(srow)
            urow[1] = srow[1]
            urow[2] = srow[2]
            urow[3] = srow[3]
            urow[4] = srow[4]
            urow[5] = srow[5]
            urow[6] = srow[6]
            urow[7] = srow[7]
            urow[8] = srow[8]
            urow[9] = srow[9]
            urow[10] = srow[10]
            urow[11] = srow[11]
            urow[12] = srow[12]
            urow[13] = srow[13]
            urow[14] = srow[14]
            urow[15] = srow[15]
            urow[16] = srow[16]
            urow[17] = srow[17]
            urow[18] = srow[18]
            urow[19] = srow[19]
            urow[20] = srow[20]
            ucursor.updateRow(urow)
            n+=1
    print(urow)

The problem is that it is extremely slow, so I would like to speed it up by stopping the search cursor loop after it finds the result. The values are all unique so after it returns a result, there is no need for the loop to keep on going. That was my thought on speeding it up, but of course I am open to other suggestions, or approaches that could help me achieve my goal.


texts: Changing srs-definition string in geopackage with ogr2ogr

texts:
When creating a geopackage with ogr2ogr in WGS84 (EPSG:4326), the srs-definition string cannot be interpreted correctly when loading the data into ArcGIS. ArcGIS does not project the data corretly into my local SRS. I found, that ArcGIS expects another string for the SRS. So I am trying to set the new string with ogr2ogr and using the -a_srs -Option, as suggested here:

arcgis-geopackage-coordinate-system

Here is my command:

ogr2ogr -f GPKG test2.gpkg -a_srs "GEOGCS[^"GCS_WGS_1984^",DATUM[^"D_WGS_1984^",SPHEROID[^"WGS_1984^",6378137,298.257223563]],PRIMEM[^"Greenwich^",0],UNIT[^"Degree^",0.017453292519943295]]" test.gpkg

Inside the geopackage there is one line-layer, called "line_4326_gp".
What am I doing wrong?

I retreive the string with this command:

gdalsrsinfo -o wkt_esri epsg:4326


texts:
I have an ASCII Grid file + .prj file and can visualize this in qGIS.  I can colorize it by going to properties and using a Pseudocolor Ramp.  However, I want to export this now as a TIFF / Raster.  Any ideas on how I can properly export this?  Right clicking the layer, there is no "save-as" option :/

I assume I need to somehow apply the color ramp to the ascii grid first, instead of just visualizing the color ramp.

Thanks


texts:
I need to use a waste dump .dxf on a map but unfortunately I'm not able to match the projection with the others I have. I'm working with Google hybrid using CRS: EPSG:32632 WGS84 / UTM zone 32N and this is the area I'm working on:

And I'm trying to use this .dxf which is projected far away from the ara of interest:

I think there's something wrong with the projection I'm using, but even changing the CRS of the layer doesn't change the result.

Here's the dxf: https://infinit.io/_/3cygVkS


texts: Pixel values in a mosaic with TIFF images in Photoscan?

texts:
I´m trying to understand how works the creation of a mosaic in photoscan or in other photogrammetry software.

When I use various TIFF images with information in their pixels (in my case temperatures) to make a mosaic, how the software works when stitch them?

Does it create an average of the value of same pixel in different images?

Obviously the temperature in same pixel in different images can´t be the same, because of that I need to know what the software do when stitch the images


texts:
I want to create an application platform based on PostGIS data. I read the documentation where I saw the prerequisite to install GEOS, GDAL and Proj4. The main problem is I am using windows platform.
Can anybody refer me where I can get full documentation on how to install PostGIS, GDAL, GEOS, Proj4 on Windows with all their libraries?


texts:
How can I identify this coordinate system? When I load the .shp with QGIS I get the correct linear mapping. However when I try to convert to GeoJSON and plot with something like Leaflet.js I get nothing. These coordinates exist roughly around: 35.161066,-114.576794

[(955999.8911220022, 26397850.011872), (955998.9291220009, 26397799.942871995), (955998.9461220019, 26397798.020871997)]


texts: Invalid geometries in QGIS 3 seems to be more prevalent than QGIS 2

texts:
I have a large polygon vector data set containing approximately 10,000 rows of data. I have used this data to run specific tasks in the QGIS 2 series (many of the Geoprocessing tools) with no issue of invalid geometries.

In QGIS 3 it seems to be picking up many more geometry validity errors even though I am using the exact same data set. I have tried multiple times in both QGIS 3 and 2 on different samples and QGIS 3 can sometimes encounter errors when 2 does not.

I know you can set the parameters to ignore features with invalid geometries but this is leaving me with a lot of errors and blanks in my analysis. Another thing to note is that this is the case in all updates to QGIS 3 up to 3.2.1.

Does anyone know why this might be happening in QGIS 3 and if there is an easy fix?


texts: Obtaining specific part of string from field in QGIS attribute table?

texts:
I'm having a bit of trouble here with this situation.

I want "NumeroCamp" only to have the numbers are showing in "Text".

Is there any function in field calculator that could help me with that?


texts: openlayers extent intersects is not working properly

texts:
I am using

ol.extent.intersects(extent1, extent2)

to check for intersection between layers and polygon I drew, but it doesn't get me any intersection which is weird because I draw the polygon creating intersection with the layer

can you think of any reason of this please ?


texts: OpenLayers: programmatically fire ModifyFeature event

texts:
Does anyone know how to programmatically simulate a vector feature click that puts this feature from normal mode:

to modifyfeature mode:

I only have one feature inside the vectors layers and i have a edit button. I have everything working and i'm just trying to save one unnecessary feature click.


texts: Transfering features class from ArcSDE to File Geodatabase using ArcPy?

texts:
I would like to know if it's possible to use python to transfer features class from SDE into a GDB and then transfer it back to SDE without losing SDE parameters such topology, versionning, relationship.

The reason of this loop is that I would like to use FME to edit the features class, but I dont have the license to read/write into SDE directly so I would like to do the updates inside a GDB.

Can somebody tell if it's possible, and if yes, expose the arcpy functions that I shall use?

I have this so far

Before running the FME translation:

#import FeatureCLass from SDE into GDB used as input for FME

import arcpy

SOURCE = "C:/Users/DemMa18/AppData/Roaming/ESRI/Desktop10.3/ArcCatalog/GEOP011@demma18.sde/DEMMA18.TestFME"
DEST = "C:/temp/arcgis/FME/TestFME_IN.gdb"

#delete all FC in GDB
arcpy.env.workspace = DEST
fcDestList = arcpy.ListFeatureClasses()
for fc in fcDestList:
    arcpy.Delete_management(fc)

#copy SDE FC into GDB
arcpy.env.workspace = SOURCE
fcSourceList = arcpy.ListFeatureClasses()

for fc in fcSourceList:
    arcpy.Copy_management(fc, DEST+"/"+fc)

After FME translation:

#export FeatureCLass from GDB used as output by FME to SDE

import arcpy

DEST = "C:/Users/DemMa18/AppData/Roaming/ESRI/Desktop10.3/ArcCatalog/GEOP011@demma18.sde/DEMMA18.TestFME"
SOURCE = "C:/temp/arcgis/FME/TestFME_OUT.gdb"

#Delete features in SDE FC
arcpy.env.workspace = DEST
fcDestList = arcpy.ListFeatureClasses()
for fc in fcDestList:
    arcpy.TruncateTable_management(fc)

#Append GDB FC into SDE FC
arcpy.env.workspace = SOURCE
fcSourceList = arcpy.ListFeatureClasses()

for fc in fcSourceList:
    arcpy.Append_management(fc, DEST+"/"+fc)


texts:
I wanted to test QGIS 3.18 Parametrizable SVG functionality. More precisely I wanted to test text replacements (as described by some posters here)
So far I'm trying to make some very very simple .svg using Inkscape.
My main issue is that even with these simple SVGs the text is not visible, not well rendered or at least badly positioned. Simple SVGs such as a single letter ou a letter at the center of a circle are always inaccurate.
How can I troubleshoot this behaviour ?
Is there any path to follow to maximize the QGIS compatibility of the generated SVGs?
Is Inkscape a tool to avoid?
I really would like to avoid creating SVG files by hand or starting from scratch.


texts:
I would like to script QGIS to import a style library as part as a configuration for a new user.  By Python commands, how do I instruct QGIS to  do what a QGIS user can do with the "Import Items" dialog box that is within the "Style Manager".  Specifically, I would like to instruct QGIS to: Import from "File" (not "URL"); check to "Add to favorites"; blank for "Additional tag(s) and "Select All".

When I look at the QGIS Python API I have found classes for working map layers or for prompting the Style Manager Dialog.  From my understanding the map layers requires a layer to exist and isn't suitable for a blank workspace, whilst, the Dialog option doesn't allow one to set all of the options.

Some pseudo code could be like:

newStyleCollection = QgsStyleManager()
newStyleCollection = newStyleCollection.importXML('c:\\SymbolSet.xml')
newStyleCollection.AddToFavourites('True')
newStyleCollection.AdditionalTags(None)
newStyleCollection.sync()

I am basically after the libraries for working with the QGIS style manager.


texts: Get zoomed-in images with the EarthExplorer API

texts:
I've been using the landsatxplore library that provides a Python API for EarthExplorer. However, when I've been searching for scenes using code like this:

scenes = api.search(
    dataset='LANDSAT_ETM_C1',
    latitude=40.906,
    longitude=-73.810,
    start_date='2017-01-01',
    end_date='2019-01-01',
    max_cloud_cover=0.5,
)

the only scenes that have been returning are ones that encompass half of the northwest USA, like this:

When in reality I'm looking for images that show individual houses and roads, like this (direct from the EarthExplorer interactive map):

Its obviously technically possible to get images like the above, because EarthExplorer possesses these qualities of images, but none of the scenes that have been returning with my code have them. Is there a way, through the EarthExplorer API and landsatxplore, to get images that are zoomed in like this? Barring that, is there an API of some sort that, when given the coordinates of a location, will provide a similar image?

EDIT: To clarify, I'm looking for services that can provide images with the same or slightly lesser zoom levels, as shown in the second image, of the entire world, or close to it, not just the US.


texts: Spatial join for closest features and using a sum merge rule in QGIS

texts:
I am trying to perform a spatial join based on nearest neighbor ("CLOSEST") spatial relationship in QGIS, and by summing the results of the join where multiple matches occur.

In ArcMap there is a way to do this using the spatial join tool, selecting the "CLOSEST" Match Option, and specifying a "SUM" merge rule for the joined attributes.

I would like to find a QGIS equivalent of this operation. I have already found a plugin which can do a "nearest neighbor join" called NNJoin. I am reasonably satisfied with the results of this plugin. Please note my input shapes are all MultiLineStrings (both target and join feature).

Now I would like to join this resulting table from NNJoin back to my target feature, but it has to sum things for multiple matches. I started fiddling around with the Virtual Layer option in QGIS in order to come up with some SQL code that can do this for me, but I have the impression that the Virtual Layer options are still defective in QGIS (see Virtual layer functionality does not seem to support memory layers).

Is there any reasonably straighforward way to achieve this using QGIS?


texts: How to get the nearest point on a linestring to a given point?

texts:
I've been using PostGIS for a long time now but never had to use the LINESTRING geometry ...! :)

Here's what I would like to do : I have a table of linestrings (representing streets of a given city, SRID 3395)  and I would like to find the nearest linestrings to a given point (GPS position, SRID 4326).

The solution I found is to select all the linestrings within my point using the expand() method and determine the distance between each linestring and my point using the ST_Distance() method.

Here's the SQL :

SELECT myLineId, myLineName, ST_Distance(ST_Transform(GeomFromText('POINT(LON LAT)',4326),3395),myLineGeom) AS myLineDistance
FROM myLines
WHERE myLineGeom && expand(ST_Transform(GeomFromText('POINT(LON LAT)',4326),3395), 100)
ORDER BY myLineDistance;

The results I get look OK but I have the feeling something is wrong in my implementation.

1) Do you guys think the expand() can  get all the concerned linestrings?

2) Do you guys think the ST_Distance() is the right method to use? I guess I'm doing it wrong since the distance I would like to get is the smallest distance between the point and my line and not the distance between the point and one of the points of the linestring.

Illustration:


texts: GeoTools JTS: Snap a Point to a polygon (equivalent to PostGis ST_ClosestPoint)

texts:
I'm using GeoTools JTS libraries in Java and I want to find the nearest point inside polygon shape to some single target location. This is exactly what in PostGis ST_ClosestPoint does.

I found this article discussing "Snap a Point to a Line".

So trivial solution will be:

-Dump polygon to lines
-Snap a Point for each above Line
-Calculate the distance between the point and "the snapped" point
-Take the one with minimum distance

How can I achieve my requirement in GeoTools in more "elegant" way?

EDIT:

Actually I'm not sure if the link to snap a point doing exactly what ST_Closest because I don't understand how it's actually extracting the closest point...


texts: QGIS LYON desktop cannot be opened on windows 7

texts:
I have installed QGIS lyon from

http://qgis.org/downloads/QGIS-OSGeo4W-2.12.3-1-Setup-x86_64.exe

on win7.

I did not install the default data set.

On my desktop, there are some icons of QGIS.

When I clicked "QGIS Desktop 2.12.3" icon, I got a window telling me that it is setting up GUI and checking provider plugins.

But, after that, it disappeared. Nothing opened on my desktop.

Did I miss something ?


texts: QGIS export 3d view / vector layers of 3d model to cad/blender friendly format

texts: I am loading vector layers in .gml format that create 3d models on a map. I'd like to save that model out of those layer using either 3d rendered map or other method. Is it possible, and If so, how?


texts: Grayscale Bing Basemap - Is it possible in ArcMap?

texts:
Is it possible to achieve a grayscale look for the Bing basemap in ArcMap 10?  I noticed that its possible with the maps hosted by ESRI since we can access the server that is hosting them.  I haven't been able to find anything on Bing.  Even using the OpenStreet Maps would be a great alternative, but again there doesn't seem to be an open server to use the trick posted on ESRI's mapping blog.

Edit

It looks like this is currently not possible, so I'm accepting the "its not possible answer".  If this does become possible in the future please submit an answer!


texts:
I want to create a tool layer in QGIS 2.4  that will select any vector features that intersect with that point. In my script I tried the following function:

def selectFeature(self, point, button):

       pntGeom = QgsGeometry.fromPoint(point)
       pntBuff = pntGeom.buffer( (self.canvas.mapUnitsPerPixel() * 2),0)
       rect = pntBuff.boundingBox()
       cLayer = self.canvas.currentLayer()
       selectList = []
       if cLayer:
               provider = cLayer.dataProvider()
               feat = QgsFeature()
               provider.select([],rect)
               while provider.nextFeature(feat):
                       if feat.geometry().intersects(pntGeom):
                               selectList.append(feat.id())
               cLayer.setSelectedFeatures(selectList)

But I have the following error:

  AttributeError: 'QgsVectorDataProvider' object has no attribute 'select'
  in line "provider.select([],rect)"

What am I doing wrong?


texts: Calculating latitude and longitude of points using QGIS

texts:
I'm trying to figure out how to calculate the lat/long of points using the QGIS Field Calculator.
Is there an easy way to do this?
I saw a previous related question (Using QGIS API and Python, to return latitude and longitude of point?), but it deals with a more advanced approach.
In ArcGIS Desktop this is similar to the Calculate Geometry function; I just can't find where do this in QGIS.


texts: Checking "connectedness'' of line shapefile in ArcMap?

texts:
I inherited a line shapefile in ArcMap that is used to represent a road network.  The problem is that the road network is pretty large and it is impossible for me to:

Tell if the end of each line (road segments) is "connected" to either an edge, vertex, or the end of another line (road segment); and
Tell if the entire road network is "connected" - i.e., can I trace a route from any point on the road network to any other point on the road network?

Is there a way to accomplish this in ArcMap?  That is, is there a function that will allow me to highlight line segments that do no connect to other line segments, or something similar that would allow me to answer the two questions above?


texts: Adding message bar to ArcMap using ArcObjects?

texts:
Is it possible to create a message bar popup(present in QGIS) in ArcMap, instead of shooting users with MessageBoxes?

I am not expecting any code examples. Guide me on the approach if this is possible.


texts: Optimizing ST_Within query to count lightning occurrences inside country

texts:
I just started learning PostGIS, and right now i'm trying to count the
number of lightning strikes in my country (Brazil in this case).

I have the following tables:

earth_feed --> This holds the geometry(Point, 4326), lon, lat and some other info about a specific lightning occurrence.

brasil_shape --> This table has a single row geometry(MultiPolygon, 4326) which is the outfit of Brazil.

My query is for 14-January (the whole table is for Jan/2014, it's a test table)

SELECT COUNT(earth_feed.geom) FROM
earth_feed as earth,
brasil_shape as brasil
WHERE
EXTRACT('day' FROM date_occur) = 14 AND
ST_WITHIN(earth.geom, brasil.geom)

This query takes about 5 minutes to complete, and for this day we had about 200k lightning occurrences. Also, the table indexes seems to be ok:

CREATE INDEX idx_light_geom
  ON earth_feed
  USING gist
  (geom);

So i have 2 questions:

1 - Is this a 'normal' execution time? Brazil has continental dimensions and a lot of lightning occurrences.

2 - It seems that the shapefile for the country outline is being 'pulled' at each comparison. Here's part of the EXPLAIN output:

->  Nested Loop  (cost=0.00..9.58 rows=1 width=32)"
     ...
     ->  Seq Scan on brasil_shape brasil  (cost=0.00..1.01 rows=1 width=32)

This affirmation is correct?

3 - The country outline is a MultiPolygon, so what happens if i change this shapefile to one having a smaller grid (more polygons)? Will the index be more effective? (I didn't tried this yet).

Here's the full EXPLAIN output (i changed the names to make the post easiear to read)

"Aggregate  (cost=9.59..9.60 rows=1 width=32)"
"  ->  Nested Loop  (cost=0.00..9.58 rows=1 width=32)"
"        ->  Seq Scan on shape_brasil brasil  (cost=0.00..1.01 rows=1 width=32)"
"        ->  Index Scan using idx_raios_diario_jan_2014_geom on raios_diario_jan_2014 earth  (cost=0.00..8.56 rows=1 width=32)"
"              Index Cond: (earthnetworks_geometry_raios && brasil.geom)"
"              Filter: ((date_part('day'::text, earthnetworks_dt_horario) = 14::double precision) AND _st_contains(brasil.geom, earthnetworks_geometry_raios))"


texts: Converting ArcMap shapefile with labels to KML (labels retained)

texts:
I am having problem i converting a shapefile with label to Google Earth KML using ArcGIS Desktop 10.1. I would like to convert a shapefile with label into KML wherein the label is retained in KML and readable. I have tried to convert the shapefile but the label in kml is unreadable.

.


texts: Can somebody help me in making 3-dimensional map from 2-dimensional??some software reference or some tool which can do this. I know the height of the building and other aspects....