texts: Issue with arcpy.MakeFeatureLayer_management in tool

texts:
I have created a python toolbox that creates separate layers based on unique field values within a feature layer.  For some reason arcpy.MakeFeatureLayer_management will create the layer but will not add it to the TOC even after using RefreshTOC().  I used the same tool in the python interpreter and it adds the created layer to the TOC.

The layer doesn't need to be written to disk or saved permanently; I wanted to leave the saving up to the user.  Any thoughts on why there is a difference and how to get the layer added to the TOC?

EDIT

    ws = arcpy.env.workspace
    mxd = arcpy.mapping.MapDocument("CURRENT")
    data_frame = arcpy.mapping.ListDataFrames(mxd)[0]
    layer = parameters[0].valueAsText #GPFeature Layer
    field = parameters[1].valueAsText #field name

    seach_layer = arcpy.mapping.Layer(layer)
    lfields = arcpy.ListFields(seach_layer.dataSource)

    rows = arcpy.SearchCursor(layer, fields=field)
    #get they field type to use in query
    ftype = None
    for f in lfields:
        if f.name == field:
            ftype = f.type
    # get values for unique list
    value_list = []
    for row in rows:
        value = row.getValue(field)
        value_list.append(value)

    unique = set(value_list)
    len_list = len(unique)

    try:
        for name in unique:
            if ftype in ['Float','Double','SmallInteger']:
                where_clause = "\"%s\"=%s" %(field, name)
            else:
                where_clause = "\"%s\"='%s'" %(field, name)

            temp = r"%s_temp_%s" %(layer, field)

            # make a feature layer to an "outfile" that's in memory
            # add to TOC and leave up to user to save if they want
            arcpy.MakeFeatureLayer_management(layer, temp,
                                              where_clause)

    except Exception as e:
        arcpy.AddError(e)


texts:
I have an MXD that I would like to publish to an ArcGIS server using C# ArcObjects (preferably as an MSD). I have been unable to find anything regarding this. Essentially I would like to implement the following ArcPy code in ArcObjects:

arcpy.mapping.ConvertToMSD()
arcpy.mapping.PublishMSDToServer()

Any help would be appreciated! Thanks!


texts:
I have a layer that gets created with a couple of attributes, one of which is NAME. Then I create a LabelLayer and add a feature layer like so:

layerLabel.addFeatureLayer(layerDisplay, rendererLbl, "{NAME}", { pointPriorities: 'CenterCenter' });

And this works fine. The NAME is displayed where I'd expect it to be displayed. But, I also want to add more attributes, so I have some code that adds an attribute called VALUE and then, on the update-end event of layerDisplay, I add my label layer like so:

layerLabel.addFeatureLayer(layerDisplay,  rendererLbl, "{NAME} {VALUE}", { pointPriorities: 'CenterCenter' });

But it will display as, for example: Alaska {VALUE}. Why won't it see the VALUE attribute? If I add a break point before the label layer is added, I can do this:

console.log(layerDisplay.graphics[0].attributes);

And I can see that both the NAME and the VALUE attributes are there.

A little hacking around, and I could (but would rather not) do this before adding my feature layer:

layerLabel._buildLabelText = function (b, d, c, e) {
    return d.NAME + " " + d.VALUE;
}

Overwriting the existing internal _buildLabelText function. Obviously, I'd prefer not to have to do this, but it is instructive that the object passed as d has all the properties (including NAME and VALUE) that I want, so I'm still confused as to why the regular function was unable to match them up.


texts: Problem with nodata when exporting only a small portion of a VRT with multiple image overlaping

texts:
I have a bunch of aerial photographs (>70 covering over 10000km2) that are tiled from which I want to build a VRT and ultimately export multiple (~2000) tiny sections (~300m x 300m) distributed around the whole territory using gdal_translate.  That should be fairly easy, however I have some tiles that are represented twice (for 2 differents years) where the difference is in the NODATA.  Here is a simple example:

Image 1:

Image 2:

When I build the VRT using:

"C:\Mrnmicro\Applic\gdal\1.11.1\gdalbuildvrt.exe" -srcnodata "0" -overwrite  "C:\Users\ferba1\Desktop\jeter\testvrt\mosai.vrt" "C:\Users\ferba1\Desktop\jeter\testvrt\*.jp2"

I get a nice result:

However, if I try to export the two squares using gdal_translate:

C:\Mrnmicro\Applic\gdal\1.11.1\gdal_translate.exe" -projwin 273000 5246300 273300 5246000 -of "GTiff" "C:\Users\ferba1\Desktop\jeter\testvrt\mosai.vrt" "C:\Users\ferba1\Desktop\jeter\testvrt\nord.tif
C:\Mrnmicro\Applic\gdal\1.11.1\gdal_translate.exe" -projwin 278000 5238300 278400 5238000 -of "GTiff" "C:\Users\ferba1\Desktop\jeter\testvrt\mosai.vrt" "C:\Users\ferba1\Desktop\jeter\testvrt\sud.tif

The red square is all nodata while the blue square is ok.  Is there a way that gdal_translate on a vrt with no data like that better deal with missing data when I have overlapping images?


texts: Counting points on a line or polygon using join: number of points are calculated incorrectly

texts:
I am attempting to count the number of stream crossings within approx. 100m segments of a trail. I used this article to get quite close with join:

http://support.esri.com/cn/knowledgebase/techarticles/detail/30779

Unfortunately, this method fails to count a point once and awhile. For example, it counts 2 points on the line when in fact there are 3. I have tried creating a buffer on the points and the result is identical, so I don't think it is a proximity issue. I have also tried creating buffer on the line.


texts: Clip feature in QGIS doesn't return desired result

texts: I wanted to georeference a JPG map to get the shapefile of a city. I first uploaded the shapefile file of the entire state which had all the cities. then I georeferenced an existing JPG map of the state which had all the boundaries of all the cities in the states. Afterward, I created a new empty shapefile and edited it with polygons so I could get the boundary of the city I am interested in. The last thing I had to do was to clip the new file but I when clipped it I didn't get the desired results, as illustrated in the diagrams showing entire clipping procedure.


texts: Create prj4 file for NAD27 Mexico from prj file for import to QGIS

texts:
I would like to be able to add the NAD27 MEXICO projection to QGIS, MapInfo, etc. Global Mapper has this projection, and the prj file is:

PROJCS["NAD_1927_UTM_Zone_13N",GEOGCS["GCS_North_American_1927",DATUM["D_NAD27_MEXICO",SPHEROID["Clarke_1866",6378206.4,294.9786982]],PRIMEM["Greenwich",0],UNIT["Degree",0.017453292519943295]],PROJECTION["Transverse_Mercator"],PARAMETER["latitude_of_origin",0],PARAMETER["central_meridian",-105],PARAMETER["scale_factor",0.9996],PARAMETER["false_easting",500000],PARAMETER["false_northing",0],UNIT["Meter",1]]

Can anyone help me make this into an accurate, working prj4 file so I can add it as a custom projection in QGIS?


texts: Finding census block with highest pop. value within census tract using QGIS?

texts:
In QGIS, I'm trying to plot a centroid of only the census block with the highest population within each census tract. So there should only be one point per tract- that of the block with the greatest population. The two shapefiles that are relevant here are: 1. Census Tracts, which has the ID of each tract in the attribute table, and 2. A population density grid at the census block level, which has the pop. count per block, the tract ID, and the blocks ID in the attribute table. See screen shot. Thick lines are tracts, thin lines are blocks.

The end result is I want to do a nearest hub analysis based on the block with the highest population within any given tract.

I'm new to GIS, so you may have to spell some things out for me. I'm only about a month in.


texts: Obtain percentiles from an image in Google Earth Engine

texts:
Suppose I have a 1-band image that represents a species' predicted density and distribution. I'd like to be able to spatially define regions based on density. For instance I'd like to extract pixels that represent the top 10% of of the population. That is, identify the largest valued pixels that when summed are 10% of the total sum of all pixels. Or top 20%, top 30% and so on.

I don't know if this is the best way but in my mind I see the workflow as the following possibility: 1) convert image to array, 2) sort those values from min to max and calculate a cumulative sum, 4) convert back to an image and divide by the maximum value...now any pixel value greater than .9 represents a pixel in the top 10%.

This is the script I've been playing with.

https://code.earthengine.google.com/23ea182e407ed56f8a822ca85dd80d31


texts: Field-based hyperlinks in an ArcMap map package (.mpk)?

texts:
I have a point shapefile where each point has multiple hyperlink fields in its attribute table as per How to add multiple relative hyperlinks to single point in ArcGIS for Desktop with Basic level license? These hyperlinks are to photos stored locally.

Ideally I want to bundle up the map (shapefiles and all) and the photos and export them as a map package (*.mpk) so that the receiver can open it up and click on the hyperlinks in the shapefiles.

Is there any way of doing this? I know the map package export function has the option to include additional files but I am not sure how to get this to work for me.
Attachments are not an option, as my licence is Basic, and geodatabases or dynamic hyperlinks are not options as I need each point feature to link to multiple photos.


texts: Compute forest loss for polygons in Google Earth Engine

texts:
My query is in two parts:
Part 1:
I am trying to calculate forest loss per polygon. I uploaded coordinate points to GEE and then I drew a polygon around the points. I created a buffer of 250m to the polygons. After that I then computed forest loss per the buffer feature collection I created but when I run the code it just returns an overall figure rather than forest loss per polygon.
// Make a list of Features.
var features = [
  ee.Feature(Diamante),
  ee.Feature(Romero),
  ee.Feature(Shintuya),
  ee.Feature(MLC)
];

// Create a FeatureCollection from the list and print it.
var fromList = ee.FeatureCollection(features);
print(fromList);

// Create a geodesic polygon.
var polygon = fromList.geometry()

// Compute a buffer of the polygon.
var buffer = polygon.buffer(250);

// Compute the centroid of the polygon.
var centroid = polygon.centroid();
Map.addLayer(buffer, {}, 'buffer');
Map.addLayer(centroid, {}, 'centroid');

// Get the forest loss image.
var lossImage = gfc.select(['loss']);

// Sum the values of forest loss pixels in the Congo Republic.
var stats = lossImage.reduceRegion({
  reducer: ee.Reducer.sum(),
  geometry: buffer,
  scale: 30,
  maxPixels: 1e9
});
print('pixels representing loss: ', stats.get('loss'));

// Set the scale for our calculations to the scale of the Hansen dataset
// which is 30m
var scale = gfc.projection().nominalScale();

// Create a variable for the original tree cover in 2000
var treeCover = gfc.select(['treecover2000']);

// Convert the tree cover layer because the treeCover by default is in
// hundreds of hectares, but the loss and gain layers are just in hectares!
treeCover = treeCover.divide(100);

// Create a variable for forest loss
var loss = gfc.select(['loss']);

// Create a variable for forest gain
var gain = gfc.select(['gain']);

// The units of the variables are numbers of pixels
// Here we are converting the pixels into actual area
// Dividing by 1 000 000 so that the final result is in km2
var areaCover = treeCover.multiply(ee.Image.pixelArea())
                .divide(1000000).select([0],["areacover"]);

var areaLoss = loss.gt(0).multiply(ee.Image.pixelArea()).multiply(treeCover)
              .divide(1000000).select([0],["arealoss"]);

var areaGain = gain.gt(0).multiply(ee.Image.pixelArea()).multiply(treeCover)
              .divide(1000000).select([0],["areagain"]);

// Sum the values of loss pixels.
var statsLoss = areaLoss.reduceRegions({
  reducer: ee.Reducer.sum(),
  collection: buffer,
  scale: scale
});
// Sum the values of gain pixels.
var statsGain = areaGain.reduceRegions({
  reducer: ee.Reducer.sum(),
  collection: buffer,
  scale: scale
});
Export.table.toDrive({
  collection: statsLoss,
  description: 'NP_forest_loss'});

Export.table.toDrive({
  collection: statsGain,
  description: 'NP_forest_gain'});

Part 2: I tried to download the polygons as shapefiles but when I tried it didn't work, is there a piece of code to take the whole feature collection and then download individual shapefiles rather than just one shapefile with all locations?
Here is a link to my code:
https://code.earthengine.google.com/?scriptPath=users%2Ftravelmark85%2Fg_engine_course%3Aforest_loss_activity


texts: Sentinel-2 Tiles - Find a specific point in the image

texts:
I am new to the area of Remote Sensing. I am working with the Sentinel-2 data from AWS. I've downloaded this specific tile and created an RGB GeoTiff. Now, I want to crop the sea in the bottom of the image. First, I get a specific longitude and latitude (here: 42.364456, 0.183105). Since gdalinfo tells me that the coordinate system is "WGS 84 / UTM zone 30N" I converted the longitude and longitude to UTM which results in the zone 31T. But this does not work with the coordinates of the image (it is tile 30T, I already know that the Sentinel tiles are overlapping) How can I solve this problem? Longitude 42.364456 and latitude 0.183105 result in UTM 31T 268039 4694087.

In general: How can I find the exact position (longitude and latitude given) in a tile given like the example here?

Here is the output of gdalinfo of the example tile:

Driver: GTiff/GeoTIFF
Files: 30_T_YN_2015_12_3.tif
Size is 10980, 10980
Coordinate System is:
PROJCS["WGS 84 / UTM zone 30N",
    GEOGCS["WGS 84",
        DATUM["WGS_1984",
            SPHEROID["WGS 84",6378137,298.257223563,
                AUTHORITY["EPSG","7030"]],
            AUTHORITY["EPSG","6326"]],
        PRIMEM["Greenwich",0,
            AUTHORITY["EPSG","8901"]],
        UNIT["degree",0.0174532925199433,
            AUTHORITY["EPSG","9122"]],
        AUTHORITY["EPSG","4326"]],
    PROJECTION["Transverse_Mercator"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",-3],
    PARAMETER["scale_factor",0.9996],
    PARAMETER["false_easting",500000],
    PARAMETER["false_northing",0],
    UNIT["metre",1,
        AUTHORITY["EPSG","9001"]],
    AXIS["Easting",EAST],
    AXIS["Northing",NORTH],
    AUTHORITY["EPSG","32630"]]
Origin = (699960.000000000000000,4800000.000000000000000)
Pixel Size = (10.000000000000000,-10.000000000000000)
Metadata:
  AREA_OR_POINT=Area
Image Structure Metadata:
  COMPRESSION=JPEG
  INTERLEAVE=PIXEL
Corner Coordinates:
Upper Left  (  699960.000, 4800000.000)
Lower Left  (  699960.000, 4690200.000)
Upper Right (  809760.000, 4800000.000)
Lower Right (  809760.000, 4690200.000)
Center      (  754860.000, 4745100.000)
Band 1 Block=10980x16 Type=Byte, ColorInterp=Red
Band 2 Block=10980x16 Type=Byte, ColorInterp=Green
Band 3 Block=10980x16 Type=Byte, ColorInterp=Blue

Tile:


texts: Ideal minimal area for watershed automatic delineation using DEM SRTM

texts: I need to do a bunch of analysis on a city using for that, watersheds . As smaller the watersheds are, better are quality of my work. Due to the vertical accuracy of SRTM DEM, I know there are some limitation of area. For example, as smaller the areas are, less diluted the errs are. I was wondering if anybody could shed a light on how small my watershed areas could be and still be Trustworthy. Indication of any paper or article would also help.


texts: Problem getting SAGA activated in Ubuntu QGIS 2.18.2

texts:
By some reason SAGA does not appear into Processing Toolbox list although I have checked it to be active. See screen capture. What could be the problem?

QGIS vesion is 2.18.2 and OS is Ubuntu 16.04 (Xenial).


texts: ArcMap 10.2 Shapefile Labels Don't Work With Pigiarniq Font

texts:
I've recently asked the question to the ESRI help team and apparently it's an unreported bug. I've also had a few suggestions from some people on GeoNet but none have worked. So here it goes:

The problem I am having is that my "XY Event Source" labels do not work when I convert them to shapefiles.

When I click "Label Features", my points in the shapefiles display a bunch of question marks (?).

I am using the Pigiarniq font which I have downloaded on my computer.

Even when I change the "Text Symbol" from the "Layer Properties" to Pigiarniq font instead of standard Arial, it comes up as question marks (?).


texts:
I procure data from DigitalGlobe (Quickbird Satellite data) for different periods hence data having some different band combination. I have attached an example image to show the problem.

I am looking for your guidance.


texts: Pre-compiled GIS software for ARM architectures

texts:
I would like to install some GIS software on my ARM laptop; QGis in particular would be useful. On x86 architectures I install packages from the much useful UbuntuGIS repository. However, this repository only makes available a restricted subset of packages for ARM.

Is there any other source of pre-compiled GIS software for ARM architectures? Or must I compile the stuff myself?


texts: Can't display XY data in MXD on ArcGIS Desktop- onEnable error

texts:
When using Display XY data function on imported excel into MXD I get this error:

and if I click on "ok" then I get this message:

and the MXD get stuck and I must shut down the ArcGIS program. I know this solution http://desktop.arcgis.com/en/arcmap/10.3/manage-data/editing-fundamentals/data-source-does-not-have-an-oid-field.htm

but I don't understand what has changed from yesterday? by routine I always managed to convert Excel file into shapefile with Display XY data function

I work on this Excel, all the field are regular- as always:


texts:
I need to do a simple merge polygon in PostGIS. The attribute table from both input shapefiles needs to be preserved. This is easily done In QGIS by using the merge polygon tool. But I don’t know how to do this in PostGIS.
I want to merge poly1 and poly2 – see picture (done in QGIS).

I have tried with the code below, but this gives way to many rows (should be around 5000, but it gives more 1 millions rows:
CREATE TABLE land_merge as
SELECT ST_Union(poly1.geom, poly2.geom) as geom
FROM poly1, poly2;


texts: First and last point of MultiLineString-objects in QGIS 3

texts:
I would like to identify the start- and endpoints within a set of polylines and return the coordinates of both points finally.
So far i have been able to iterate the features and found out that ArcGIS's export created a MulitLineString with a number of vertices (where the number corresponds to the length of the line).
Using feature.geometry().vertexAt(0) (where feature is the current feature within a for loop) returns the first point of every feature and the print-result looks like that <QgsPoint: Point (4800744.81731882505118847 2812251.2178244274109602)> but i didn't find a way to set up a "for vertex in vertices"-style iteration to get the last point of every line.
The Python2.x-style from Can vector layer get start point and end point of line using PyQGIS with
geom = feature.geometry().asPolyline()
start_point = QgsPoint(geom[0])
end_point = QgsPoint(geom[-1])

doesn't work anymore.asPolyline() seems to return an empty list.
I've tried to iterate the vertices with vertices(self) → QgsVertexIterator but cant get any results.
since i am pretty new to PyQGIS the QGIS Python API seems really confusing to me.


texts:
I want to show the coordinates of my map in gauss-krüger coordinates, but nothing works. I tried another projection, displayProjection, lonLat.transform but the answers are the same.
my source projection is EPSG:900913


texts: How to use the TopoViewer in the QGIS DB Manager?

texts: I'm just wondering what the TopoViewer in QGIS is doing. DB Manager (PostGIS) says 'Schema "schemaname" is not registered in topology.topology'. Unfortunately there is no help available and I can't find a manual.


texts: Given a shapely.geometry.linestring: How do you access single vertices using Python?


texts: Disambiguate messy place names in python (preferably on local machine)

texts:
I have list with several million place names that come from Flickr profiles.  Users provided these placenames as free text, so they look like this:

Roma, Italy
Kennesaw, USA
Saginaw, MI
Rucker, Missouri, USA
Melbourne, Australia
Madrid, Spain
live in Sarnia / work in London, Canada
Valladolid, España
Italia
West Hollywood, United States

I want to disambiguate these place names.  I am aware that there is in some cases no straightforward to this solution, but I am willing to live with some false disambiguation and with "no answer" for some of the places.  If a place name corresponds to the name of multiple cities, then I want to assign that place to the largest city that it corresponds to.

Yahoo's place finder api would be a good solution to this problem, but I would need to make too many API calls to get through my list, so I'd like a local solution (i.e., one that does not depend on a remote api).  Does anyone know of any python libraries that do this kind of thing, or any other local solutions?

(I've also asked this question on stackoverflow.)


texts: Simple Linux desktop navigation system for tracking a drone

texts:
Trying to build my own ground control station (GCS) for different drones I have.  One is a boat and the other a plane.  I have looked at the GCS's available and they are great, but I want to do something much simpler.  I really only need.want to put a red dot on a GeoTiff, stored natively on the GCS.

The TDMA radio is Ethernet based and the drones have GPSd running on their Raspberry Pi, which publish to TCP port 2947.

What packages or applications would you recommend to achieve this ?

Should I look at hacking GPSDrive or QGIS.  QGIS looks great, but has lots of bells and whistles.

Is there something simpler ?

Would prefer to avoid anything written in Java, preferring C, Java Script and Python.


texts: Google Drive and Backing-up/Syncing GIS Data

texts:
Seeing that many of us have three computers these days (e.g.. work, home and laptop) is anyone using Google Drive as their main file repository for GIS data?  I like the idea of backing everything up in the cloud and syncing it to all of my computers.

I would like to hear anyone's experiences with backing up and syncing their GIS data in Google Drive or one of the other cloud sync/backup services.  Are there any serious limitations in using/storing GIS data in the cloud?  Is anyone integrating Google Drive with their work server to access data at home?


texts: CRS different after r.reclass through QGIS different than source raster

texts:
According to the properties, my output raster is in a different projection than the input after running GRASS r.reclass through QGIS. Does anyone know if this is just a fluke and I can simply select a different CRS through the properties, or do I have to "Save as" to reproject the output raster in the correct CRS?

Using QGIS 2.14.3


texts: How to project data into both a Vertical & Horizontal Coordinate System?

texts:
I have a shapefile whose Geometry is of POINTZM type.

The XY coordinate System is NAD_1983_CORS96_StatePlane_North_Carolina_FIPS_3200_Ft_US
and the Z coordinate system is North_American_Vertical_Datum_1988

I want to transform the data, so that the data is projected in a Geographic (Lat/Long) projection, with the WGS84 horizontal datum and the EGM96 vertical datum.

I can't figure out how I can do this using the Project tool in ArcToolbox, of ArcGIS for Desktop 10.3.

How do I do this?


texts:
I need to get land cover data for the entire world. I'm using corine land cover data for Europe:

http://www.eea.europa.eu/data-and-maps/data/clc-2006-vector-data-version-2

but I have no idea how to get data for the rest of the world.


texts: Mapping mean distances to water for collared animals using ArcGIS Pro

texts:
I have shapefiles of location data for around 40 collared animals as well as a shapefile layer of water points. I want to calculate and compare the distances of individuals to water sources in the dry and wet season.
I found the "Near" analysis tool in ArcGIS Pro, which allowed me to calculate the distance between a point and its nearest neighbour.
How can I go about mapping/visualising the mean distances to water points for each individual in the wet vs. the dry season?


texts: How can I add a PostGIS query result to OpenLayers?

texts:
I'm doing a query on PostGIS that returns a lot of geoemtries as results (more than 1000 which is why I can't use Geoserver CQL) - the query WHERE CLAUSE is being changed each time by the user (its some kind of filter with UI form)

What is the best way to add this dynamic query result to map in OpenLayers?


texts: Partially processing a catalog with an extent in lidR

texts:
I am using lidR to create a DEM of a very large area (over 1 million hectares). I have already processed a portion of this area.

I would like to be able to read the entire LAScatalog so that the tiles I still need to process can access bordering tiles and use them as buffers, but don't know how to avoid processing the rest of the catalog.

I have been looking at the partial processing option as well as the catalog_intersect() function. I can make a raster of the area that I have already processed and get the extent and use catalog_intersect() to eliminate the tiles I have already processed, but that would also remove the tiles that would be needed as buffers.

Is there a way to use the ctg$processed option to select only the area not covered by a particular extent or spatial object for processing? Or is there a way to access the numbers assigned to each chunk of the catalog so that I can export the spatial object, intersect it with the processed area, and then use the ctg$processed option to select the unprocessed area?


texts: Symbology based on field names and values in ArcMap 10.5

texts:
Is there a way to apply symbology based on field names as well as the values per each field?

I have a road layer where the fields are years and each road segment has an attribute value of 1 (had work done in that year) or 0 (no work done that year).I would like to apply symbology that shows road segments with attribute value of 1 in 2018 to be blue, 2017 to be red, 2016 to be green, etc. I have fields from 2000 to 2018 so more than the 3 fields that you can bring into the 'multiple attributes' symbology option.

I know how to pick all roads that had work done and apply a symbology to that but how do I also add in the further distinction of year?


texts: Calculating lat/lons from PostGIS table which has projected coordinate system?

texts:
I have a PostGIS table in a projected coordinate system (EPSG: 28355 or MGA94 Zone 55). In QGIS (with the PostGIS manager plugin) I calculate a coordinate by performing the following SQL:

UPDATE gisprod.buildings SET easting = x(centroid(the_geom))

Now I want to add a "lat" and "lon" field to this table and calculate the values using the appropriate coordinate system, in this case that would be (EPSG:4283 or GDA94).

Is there a way I can specify in the SQL statement to calculate the x/y coord using a different coordinate system?


texts: Export selected feature/s to shapefile with open source python library

texts:
I need to be able to select one feature/record at a time within a polygon shapefile, regardless of how many shapes it contains, using open source python libraries (e.g., OGR, pyshp, fiona). Then I want to use each feature/record in a different function that requires a shapefile format. Below is what I'm attempting with pyshp:

shp = shapefile.Reader(file)  ## open and read the shapefile
shp = shp.shapes()  ## read the geometry information of each shape record
shpInPoly = []  ## empty array that will contain the features in the file
for i in shp:
   shpInPoly.append(i)

Then I'd assume I could grab the first feature in the array of appended features like this:

currentPolygon = shpInPoly[0]

And use the "currentPolygon" variable in my other function with the first record/feature in the original shapefile.

But it doesn't work. The error I get in my other function which requires shapefile inputs is "shapefile.ShapefileException: Shapefile Reader requires a shapefile or file-like object." The other alternative, per the question title (and what I'm thinking is likely best), would be to export that selected feature/record and all its content (geometry and fields) to a new shapefile, to achieve the same purpose of using that in another shapefile-requring function. The only posts I found that attempt something similar are Exporting Selected Feature to New Shapefile using ArcPy? and Splitting shapefile into many shapefiles with open source?, but I can't use ArcPy.

Any ideas?


texts: Using the OS Open Background basemap with ArcGIS JS API 4

texts:
I'm relatively new to the ArcGIS JS API, but I'm looking to use the OS Open Background as the basemap. I'm using the latest (4.0) API and I can't find a working example or decent documentation that describes the process. I can find enough for v3.16, but nothing that also applies to 4.

I'm currently trying something like this:

layerUrl = "https://tiles.arcgis.com/tiles/qHLhLQrcvEnxjtPr/arcgis/rest/services/OS_Open_Background_2/MapServer",

layer = new TileLayer(layerUrl, {
    tileServers: ["https://tiles1.arcgis.com/tiles/qHLhLQrcvEnxjtPr/arcgis/rest/services/OS_Open_Background_2/MapServer", "https://tiles2.arcgis.com/tiles/qHLhLQrcvEnxjtPr/arcgis/rest/services/OS_Open_Background_2/MapServer", "https://tiles3.arcgis.com/tiles/qHLhLQrcvEnxjtPr/arcgis/rest/services/OS_Open_Background_2/MapServer", "https://tiles4.arcgis.com/tiles/qHLhLQrcvEnxjtPr/arcgis/rest/services/OS_Open_Background_2/MapServer"]
});
map.add(layer);

...but this just causes the browser to hang indefinitely before crashing, so something isn't responding as expected or the URLs aren't correct. Any pointers?


texts:
I have a 250m resolution raster file which is needed to be clipped using a mask layer. I am trying to do it with gdalwarp through Clipper tool in QGIS but it is taking forever. I tried using R (mask function) but that too didn't help.
Are there any methods to do it with less computational burden on the system?
This is the code I was using gdalwarp -dstnodata -9999 -q -cutline "D:/Mesh/MRB/ncrb trial 2.shp" -tr 0.002083333 0.002083333 -of XYZ "D:\\Mesh\\Shangguan et al. (2016) Bedrock data\\BDTICM_M_250m_ll.tif" D:/Mesh/MRB/result.tif.
This is the datafile: ftp://164.132.26.83/data/recent/BDTICM_M_250m_ll.tif (approx 7.5 GB).
This is the gdal info data about the file:
Driver: GTiff/GeoTIFF
Files: D:\\Mesh\\Shangguan et al. (2016) Bedrock data\\BDTICM_M_250m_ll.tif
D:\\Mesh\\Shangguan et al. (2016) Bedrock data\\BDTICM_M_250m_ll.tif.aux.xml
Size is 172800, 67200
Coordinate System is:
GEOGCS["WGS 84",
DATUM["WGS_1984",
SPHEROID["WGS 84",6378137,298.257223563,
AUTHORITY["EPSG","7030"]],
AUTHORITY["EPSG","6326"]],
PRIMEM["Greenwich",0],
UNIT["degree",0.0174532925199433],
AUTHORITY["EPSG","4326"]]
Origin = (-180.000000000000000,83.999167206000010)
Pixel Size = (0.002083333000001,-0.002083333000001)
Metadata:
AREA_OR_POINT=Area
ATTRIBUTE_LABEL=BDTICM_M
ATTRIBUTE_MEASUREMENT_RESOLUTION=1
ATTRIBUTE_TITLE=Absolute depth to bedrock (in cm)
ATTRIBUTE_UNITS_OF_MEASURE=cm
CITATION_ADDRESS=shanggv@bnu.edu.cn / tom.hengl@isric.org
CITATION_ORIGINATOR=College of Global Change and Earth System Science, Beijing Normal University / ISRIC - World Soil Information
CITATION_URL=http://dx.doi.org/10.1002/2016MS000686
CONFIDENCE_INTERVAL=M
DATA_FORMAT=Byte
DEPTH=None
DOWNLOAD_FTP_URL=ftp://ftp.soilgrids.org/data/recent/BDTICM_M_250m_ll.tif
HORIZON_LOWER_DEPTH=5000.00 m
HORIZON_UPPER_DEPTH=0.00 m
KEYWORD1=bedrock
KEYWORD2=depth
NO_DATA=255
PROJECT_URL=http://soilgrids.org
PUBLICATION_DATE=2017-03-10
RANGE_DOMAIN_MAXIMUM=5000
RANGE_DOMAIN_MINIMUM=0
SERIES_NAME=SoilGrids250m
SLD_NAME=soilgrids250m:BDTICM
SLD_URL=ftp://ftp.soilgrids.org/legends/BDTICM.sld
TECHNICAL_SPECIFICATIONS_URL=https://github.com/ISRICWorldSoil/SoilGrids250m
Image Structure Metadata:
COMPRESSION=DEFLATE
INTERLEAVE=BAND
Corner Coordinates:
Upper Left (-180.0000000, 83.9991672) (180d 0' 0.00"W, 83d59'57.00"N)
Lower Left (-180.0000000, -56.0008104) (180d 0' 0.00"W, 56d 0' 2.92"S)
Upper Right ( 179.9999424, 83.9991672) (179d59'59.79"E, 83d59'57.00"N)
Lower Right ( 179.9999424, -56.0008104) (179d59'59.79"E, 56d 0' 2.92"S)
Center ( -0.0000288, 13.9991784) ( 0d 0' 0.10"W, 13d59'57.04"N)
Band 1 Block=172800x1 Type=Int32, ColorInterp=Gray
NoData Value=-32768
Overviews: 86400x33600, 43200x16800, 21600x8400, 10800x4200, 5400x2100, 2700x1050, 1350x525


texts: Drawing a square manually using coordinates

texts:
I've load points coordinates from a .csv file on QGIS. These points are the intersection of a regular grid.  I would like to select manually four point of this .csv file to draw a square polygon. Also, I would like to measure the length of the edges.
I need this square because the next step is to calculate an intersection with an irregular polygon.


texts: Creating a lightweight 1 degree x 1 degree mesh with GeoJSON?

texts:
I have produced via QGIS a 1°x1° mesh polygon and exported that in GeoJSON format. However, the file is 11-15 MB large. That is really too heavy to be used in Web applications, such as OpenLayers (which I intend to do).

The format looks like this:

{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-179,90],[-179,89],[-180,89],[-180,90],[-179,90]]]},"properties":{"left":-180,"top":90,"right":-179,"bottom":89}},
{"type":"Feature","geometry":{"type":"Polygon","coordinates":[[[-180,89],[-179,89],[-179,88],[-180,88],[-180,89]]]},"properties":{"left":-180,"top":89,"right":-179,"bottom":88}},

I wonder if there isn't a way of producing a lighter version of such a grid. Perhaps even some on-the-fly, dynamically generated layer in OpenLayers?

My aim is to lay this mesh on a global map; the user can click on a square, which sends a request to a server with the ID of the square; the server responds with a pre-calculated value for that square-ID.


texts: What does shading indicate in home range plot?

texts:
I've simulated a home range in R using the clusthr function from the adehabitatHR package. The home range is made from 50 x,y coordinates. The output is attached.

I'm trying to get clear on what the shading indicates. The shading is clearly related to the size of the polygon. Smaller polygons = darker shading, larger polygons = paler shading. But how exactly is a polygon assigned a shade? Does the shading represent a percentage of the area covered?


texts: How to label a FeatureLayer using ArcGIS Javascript 4.3

texts:
I'm struggling to determine how to label a FeatureLayer using ArcGIS Javascript 4.3.

I've found the FeatureLayer API reference here.

One of the FeatureLayer properties is labelingInfo of type LabelClass[]. The labelingInfo details says that "There is no support for labeling in 2D. Labeling is only supported in 3D SceneViews." I want to add labels to a 2D map.

I've also found the TextSymbol Class, which says it is "used to define the graphic for displaying labels on a FeatureLayer." If that's the case, I can't figure out what FeatureLayer property to include a TextSymbol.


texts: Drawing polygons from watershed boundary dataset using ArcGIS API for JavaScript?

texts:
I'm working on a web map using esri arcgis javascript api. I'd like to use the NHD Watershed Boundary Data to draw polygons based on the boundaries of the watersheds.

Does anyone have any ideas of how to start?

I've looked at the National Map where I can download data but ideally, I could pull data from an online source rather than hosting it myself where I would be responsible for any changes to the data. I'm already using the NHD Watershed Basemap as my basemap but I don't see where the watershed boundaries are available.


texts: OpenStreetMap heavily pixellated when loaded in QGIS and standalone application

texts:
I am loading an OSM street map into an application I am building. I load it as a raster layer via an XML document with the OSM GDAL data. My problem is that the map appears heavily mangled and pixellated, such that labels are virtually illegible:

I have noticed that setting the window size can sometimes make this issue less severe, but it's not much of a fix. Is this an issue with my XML file? The XML looks like this:

<GDAL_WMS>
  <Service name="TMS">
  <ServerUrl>http://tile.openstreetmap.org/${z}/${x}/${y}.png</ServerUrl>
  </Service>
  <DataWindow>
    <UpperLeftX>-20037508.3427892</UpperLeftX>
    <UpperLeftY>20037508.3427892</UpperLeftY>
    <LowerRightX>20037508.3427892</LowerRightX>
    <LowerRightY>-20037508.3427892</LowerRightY>
    <TileLevel>18</TileLevel>
    <TileCountX>1</TileCountX>
    <TileCountY>1</TileCountY>
    <YOrigin>top</YOrigin>
  </DataWindow>
  <Projection>EPSG:3857</Projection>
  <BlockSizeX>256</BlockSizeX>
  <BlockSizeY>256</BlockSizeY>
  <BandsCount>3</BandsCount>
  <Cache />
</GDAL_WMS>

As far as I know, my extents are set correctly, and reprojection doesn't seem to be relevant (I've tried it in regular Web Mercator as well as reprojecting in WGS84. It looks like this in both)

Is there something I can add to the XML document to fix this pixellation issue?


texts:
Dont kill me on CRS!

Every time I add lat longs using the same CRS as the project they are now starting to mess up and show up in Africa or Asia.... I cant seem to match up the true lat long anymore?

This never happened before I always check the project and layer projections (use a default EPSG: 3832)....

I'm sure this is a bug not a user issue?
Anyone have anything similar happen?


texts: Calculate Percent coverage using another raster, with dif resolution and orientation

texts:
I have a 0.5 meter raster that defines presence or absence (cells are either 1 or 0). I have another 5 meter raster, which I created. I would like to calculate a percent coverage of the underlying 0.5 meter raster and pull that value into the corresponding cell in the 5-meter raster. The rasters don't line up and cannot be adjusted.

One thing I now realize is that the cells are slightly less than 0.5 meters,   .46999 to be exact.    I ran block statistics, using a 10.638298 X 10.638298 rectangle (5/.469999), with snap raster set to, the slightly offset, 5-meter raster.    That's taking a while to run, but I'm guessing that it won't work quite right.   I also have a 5-meter polygon fishnet and I could also convert the .46999 meter forest cover raster into a polygon feature.   Then I could calculate the area within each 5-meter 'cell'.


texts:
I am working with a raster sub-list of raster objects (made up of sea ice concentration values in March and April over a 10 year span from NSIDC) that are projected in polar stereographic. I have cropped the raster list to a management boundary, also projected in polar stereographic, and have plotted the first item in the raster list to test if it was done correctly. Visually, it seems like it works, however I am getting weird values for the axes.

I am wondering why they are not showing up as latitude/longitude values and if there is a way to change them.

# Import shapefile
DS_boundary <- readOGR("DS_Management_Shapefile/DS_Management_area.shp")
plot(DS_boundary)

 # Change projection of shapefile to match raster list (i.e. polar stereographic)
proj4string(DS_boundary)
DS_boundary_proj <- spTransform(DS_boundary, crs('+proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +k=1 +x_0=0 +y_0=0 +a=6378273 +b=6356889.449 +units=m +no_defs'))
proj4string(DS_boundary_proj)

# Crop sub-list to the shapefile
list_mar_apr_crop <- lapply(X=list_mar_apr, FUN=crop, y=DS_boundary_proj)
proj4string(list_mar_apr_crop)

# Test crop
plot(list_mar_apr_crop$`19910301`)
plot(DS_boundary_proj, add=TRUE)


texts: Defining attributes when copying features to new memory layer with PyQGIS

texts:
I've got a question regarding the copying of features from a layer to a new layer while keeping all former attribute data. All this I want to do with PyQGIS (eventually the code will expand more).

So far I got the following code:

from PyQt4.QtCore import *

lyr = iface.activeLayer()
vl = QgsVectorLayer("Point", "temporary_points", "memory")
pr = vl.dataProvider()

pr.addAttributes([QgsField("field_1", QVariant.String),QgsField("field_2", QVariant.String),QgsField("field_3", QVariant.String),QgsField("field_4", QVariant.String)])

vl.startEditing()

outFeat = QgsFeature()
for feat in lyr.getFeatures():
    attrs = feat.attributes()
    if attrs[2] == 'pvs':
        outFeat.setGeometry(feat.geometry())
        outFeat.setAttributes(feat.attributes())
        print attrs[4]
        pr.addFeatures([outFeat])

vl.commitChanges()

QgsMapLayerRegistry.instance().addMapLayer(vl)

So far, everything works fine. However, I do not like it that I have to "hard-copy" every attribute field in the pr.addAttributes section.

What I would like to do is append all the attributes that are present in the source layer ('lyr') to the new "temporary_points" layer.
Is this possible?


texts: OpenLayers.Control.LayerSwitcher - how can I select/unselect a layer even if not available at current zoom

texts:
I love the layerswitcher feature in OpenLayers v2 (not using v3 right now).  It's great to add to the map with one line of code!  I just have one thing I'd like to change, and I'm not sure how.

When my application starts, there are defaulted visible and invisible layers. The LayerSwitcher correctly shows a checkbox if visibility is true, and unchecked if visibility is false.  Several are checked but the checkbox is "disabled" or "greyed out" because while the layer's visibility is true, I've also specified what zoom levels to turn it on/off, and if I'm outside that zoom range when I look at the LayerSwitcher, then those layers are disabled in the control.

I'd really like to be able to select/unselect layers, even if currently not displayed do to the zoom level restriction.  Maybe these could be color coded differently, but not disabled.  Perhaps there is a CSS override I could use?  Or some option I could pass in?

I suppose I could write my own layerSwitcher code, but that seems like a lot of work to change just one feature of the LayerSwitcher that comes with...

Here is an example of a layer that shows up in the LayerSwitcher but is greyed out if not at zoom levels 9, 10, or 11:

    var myLayer = new OpenLayers.Layer.WMS(
            "This layer is OFF by default",
            "http://mygeoserver.local/geoserver/gwc/service/wms",
            {
                layers : 'MyLayer',
                transparent : "false",
                styles : '',
                srs : 'EPSG:900913'
            }, {
                isBaseLayer : false,
                visibility : false,
                resolutions : map.getScales(9, 11)
            });

    map.addLayer(myLayer);


texts: How to quickly download all GMTED2010 data?

texts:
Is there a way to automatically download all the tiles in GMTED2010 data?
From:

https://topotools.cr.usgs.gov/GMTED_viewer/viewer.htm

I can only seem to be able to download by double-clicking tile by tile.


texts: Create a heatmap for VISN level data in R

texts:
I am trying to create a heatmap for Veterans Integrated Service Networks (VISN) level averages. I have found these datasets containing shapefiles.

I have the averages I want to plot based on the VISN, but I'm not sure how to go about creating this map in R.


texts: Setting fill to a variable based on data within SpatialPolygonsDataFrame?

texts:
I'm currently on the cusp of getting done what I want to get done - but having a slight problem with referring to my data from within R's sp SpatialPolygonDataFrame object. The issue is how ggplot2 deals with spatial polygons, and factors that define or are mapped to them. I have my polygons, and in the @data element, I have the relevant data I'm interested in mapping/visualising, in a SpatialPolygonDataFrame called plotData. The following code;

p = ggplot(plotData, aes(x = long, y = lat, group = id))

p + geom_polygon(fill=plotData@data$total,color='black')

yields an error about the number of elements within the aesthetic call:

Error: Aesthetics must be either length 1 or the same as the data (1651526): colour, fill

I've tried a number of different ways of accessing these data - but am not able to actually figure out how to colour it based on anything within the @data element. The difference is something to do with how ggplot2 treats polygons - (my data has only 2208 unique elements, and 2208 unique polys, not 1.6 million), but I can't seem to identify how this is put together.


texts: How to change the color of GeoJSON layer?

texts:
I want to show my road network in red color. How can I change this default color of the road network. Where should I attach style property?

var roadNetwork = new L.LayerGroup();
$.getJSON("{% url 'road_dataset' %}", function (data) {
    var geojson = L.geoJson(data, {
        onEachFeature: function (feature, layer) {
            layer.bindPopup(feature.properties.road_name);
        }
    });
    geojson.addTo(roadNetwork);
});


texts:
I am getting the following Error on using Natural neighbour. There were solutions suggested here & here. But these works fine on Windows but not on Ubuntu. I am using Ubuntu 14.04, QGIS 2.18.4 & SAGA 2.2.7.
What could be the solution for ubuntu ?


texts: Split a shapefile into smaller files on linux command line

texts: I have a shapefile which has ~350,000 entries. This is too big to handle. I'd like to split it into many smaller files, such that each file has 10,000 entries (this means there'd be ~350 shapefiles). I don't really care about what entries go to what other entries (i.e. no requirement for all with a certain attribute be in the same file, and no requirement for all entries in the same geographic area be together). The 10,000 limit is not a hard limit. So long as it's approximately 10k. I want to do this on the command line on linux (so probably ogr2ogr or similar?). How can I do this?


texts: How to get the signal when the color changes in pyQGIS?

texts:
I would like to have the signal when the color of a vector layer is changed by the user.
In the doc, there is a signal that seems to fit (http://qgis.org/api/classQgsVectorLayer.html#a403f235d3743aa59bb3aad69a72a5705) but nothing happens when it is called.

I have a method called rendererChanged:
QObject.connect(self.iface.activeLayer(), SIGNAL("rendererChanged()"), self.DoSomething)

When executed, my function DoSomething is not executed.

I have also tried :
self.iface.activeLayer().rendererChanged.connect(self.toto)

But I get this error :
AttributeError: 'QgsVectorLayer' object has no attribute 'rendererChanged'


texts: How to add more allowed values to the Headers of Geoserver

texts:
I am aware that GeoServer comes with the feature of 'Configuring HTTP Header Proxy Authentication'. However, I only understand how to customize the 'Header attribute', but the not the 'Header value'.

For example, the "sdf09rt2s: admin":

Header attribute is 'sdf09rt2s'
Header value is 'admin'

How can I create more Header value? (such as: 'user1', 'user2', etc)


texts: Imported points from CSV file (latitude, latitude) don't show up in a post code shapefile

texts:
I can produce a point vector layer by importing latitude and longitude info from a CSV file in a new project.

After saving this file as an ESRI shapefile I can import it in a Postal Code shapefile. However now, the points don't show up. The same happens if I import the CSV file immediately in the project with the postal code shapefile.

I probably do something elementary wrong. Can some one help me out?

My aim is to show locations of interest on this postalcode file.


texts: Slow plotting of spatial data in macOS RStudio's plotting device

texts:
In my macOS system RStudio, plotting spatial data (statialpolygons, points and lines) is extremely slow, making it impossible to work with.

When I try to plot the shapefile in this link (@https://drive.google.com/drive/folders/0B0XRH7FN95-RSzFRckhXNjVWYlk), it takes about 30 minutes in RStudio.
If I plot in R (outside RStudio), it takes about 20 minutes to print in an external plotting windows titled Quartz2. On the contrary, it is a couple of less than a minute in a window's machine in my lab with equivalent specs.

The best workaround I know is to use a pdf as plotting device (or png). Doing that, it takes just a couple of seconds.

shp <- readOGR("~/stof_shp", "AP_PAT_edited1")
pdf()
plot(shp)
dev.off()

But it is not a satisfactory solution. I love work with spatial data on R, but the biggest inconvenient is the difficulty for visualizing the data, and this makes it much harder.

I have talked with different people that are in my same situation, but there doesn't seem to be too much information about it when you google or search in gis.stackexchange. I have doubts on whether it affects all macOS systems, or only laptops, only some version of the OS or all of them. My system and applications are updated.

  sessionInfo()
      R version 3.4.0 (2017-04-21)
      Platform: x86_64-apple-darwin16.6.0 (64-bit)
      Running under: macOS Sierra 10.12.5

At some point, I thought it was the matter of my version of gdal or proj, but I am confident now that I have the latests set and working. I have rgdal_1.2-8 and sp_1.2-5.

> library(rgdal)
Loading required package: sp
rgdal: version: 1.2-8, (SVN revision 663)
 Geospatial Data Abstraction Library extensions to R successfully loaded
 Loaded GDAL runtime: GDAL 2.1.3, released 2017/20/01
 Path to GDAL shared files: /opt/local/share/gdal
 Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
 Path to PROJ.4 shared files: (autodetected)
 Linking to sp version: 1.2-5

I have also been messing with XQuartz and X11, getting nowhere (I wouldn't know how to get the info of which versions I have installed, but if I enter X11() on RStudio, it will try to start X11, but it will crash in an infinite loop that I have been unable to fix).

I started the question in stackoverflow, but I was told to post the question here, where it belongs better. Few months ago I made similar question here, but I did not reproduce the example data and I think that I might have been missleading in how I presented the problem (@Slow plotting SpatialPolygonDataFrames on some systems). For those reasons I am starting a new question.


texts: Openlayers 3 styling point features with icons

texts:
I'm using the following code which defines a default style, that is used to style point (and other) features in OpenLayers 3.
When I use the var image to style my features everything is Ok, and the points are rendered on my map. When I try to use var icon with the same dataset instead, nothing is displayed. In my browsers devtools I see, that the icon is loaded?

Any ideas? where the mistake lies?

//
// -------- Defining default style settings ----------
//
 var fill = new ol.style.Fill({
   color: 'rgba(255,255,255,0.4)'
 });
 var stroke = new ol.style.Stroke({
   color: '#3399CC',
   width: 1.25
 });

 var image = new ol.style.Circle({
       fill: fill,
       stroke: stroke,
       radius: 5
 });

 var icon = new ol.style.Icon({
          anchor: [0.5, 0.5],
          size: [52, 52],
          offset: [52, 0],
          opacity: 0.5,
          scale: 1.0,
          src: "./img/mmOrange.png"
 });

 var style = new ol.style.Style({
     image: image,  // this works, but when I use image: icon, nothing is displayed
     fill: fill,
     stroke: stroke
   });


texts: ApplySymbology does not recognize existing layer (ArcMap 10.6.1)

texts:
I have a script that (1) loops through a each row in a feature class; (2) creates temporary_layer_1 with the selected feature; (3) zooms to the extent of temporary_layer_1; (4) creates temporary_layer_2 with a matching field in a different feature class; (5) applies the symbology from each of the permanent (looped) layers to the corresponding temporary layers; and (6) exports the resulting map view as a PDF.
I've had success with the code on my laptop, running ArcMap 10.6.1. When I run the same code on my desktop, also 10.6.1, it fails at "arcpy.ApplySymbologyFromLayer_management," indicating that the input layer does not exist. What is puzzling is (1) that the input layer was already used successfully in the code in the previous lines, and (2) the layer types for input/output of ApplySymbologyFromLayer are identical, as one is the derivative of the other. After the processing fails, I try to run the failed line by dragging the input layer from the TOC into the Python window and it still fails. I've tried using the full file paths with no more success.
Failed to execute (Apply Symbology From Layer). Parameters are not valid. ERROR 000628: Cannot set input into parameter in_symbology_layer.
Code below:
import arcpy
import pandas as pd
import numpy as np
import os
from arcpy.sa import *
from arcpy import env

gdb = r'C:\***\***.gdb'
env.workspace = gdb
env.overwriteOutput = True

fc_2020 = 'basin_out'
fc_2010 = 'basin_x'
temp_2020 = 'basin_2020'
temp_2010 = 'basin_2010'

output_path = r'C:\**\**'

mxd = arcpy.mapping.MapDocument('Current')
df = arcpy.mapping.ListDataFrames(mxd)[0]
lyr = arcpy.mapping.ListLayers(mxd, fc_2020, df)[0]

with arcpy.da.SearchCursor(fc_2020, ['SHAPE@', 'Station']) as cursor:
    for row in cursor:
        station = str(row[1])

        arcpy.SelectLayerByAttribute_management(lyr, "NEW_SELECTION", expression)
        ext = lyr.getSelectedExtent()
        ext.XMin = ext.XMin - 100
        ext.YMin = ext.YMin - 100
        ext.XMax = ext.XMax + 100
        ext.YMax = ext.YMax + 100
        df.extent = ext
        arcpy.RefreshActiveView()

        with arcpy.da.SearchCursor(fc_2010, ['SHAPE@', 'Number']):
            EXPRESSION = '"Number" = ' + station
            arcpy.MakeFeatureLayer_management(fc_2010, temp_2010, EXPRESSION)

        # make feature with only current 2020 row
        arcpy.MakeFeatureLayer_management(fc_2020, temp_2020) # temporary layer with the selected basin

        ### change symbology for each temporary layer
        arcpy.ApplySymbologyFromLayer_management(temp_2020, fc_2020) <-- **Failure Point**
        arcpy.ApplySymbologyFromLayer_management(temp_2010, fc_2010)

        arcpy.mapping.ExportToPDF(mxd, os.path.join(output_path, station + '.pdf'), resolution = 600, image_quality = 'BEST', colorspace = 'RGB',
                                  compress_vectors = True, image_compression = 'ADAPTIVE')
        arcpy.SelectLayerByAttribute_management(lyr, "CLEAR_SELECTION")

arcpy.mapping.RemoveLayer(df, arcpy.mapping.ListLayers(mxd, temp_2020, df)[0])
arcpy.mapping.RemoveLayer(df, arcpy.mapping.ListLayers(mxd, temp_2010, df)[0])


texts: PostGIS shortest distance from point to point

texts:
I'm new to Spatial Databases and I'm trying to find the nearest point of a point.

In my application there is a google map and the user can click anywhere on it. The click executes a request to the server which should return the nearest point in the database from the point where the user clicked.

My query looks like this:

SELECT point, strength
FROM signalstrengths
ORDER BY st_distance(st_makepoint(:lng, :lat)::geography, point) DESC
LIMIT 1;

This query always returns one of two points, no matter where i click, although I have around 2k records in the database.

What is wrong with my query?


texts: Using Oracle Spatial for dynamic segmentation?

texts:
I am new to Oracle Spatial, my main objective is to conduct dynamic segmentation (DynSeg) on a selection of tables. I understand there is a rich Oracle Spatial API that has all kinds of functions and operators.

Has anyone successfully used any of the functions to implement dynamic segmentation?


texts: Alternate colors of GeoJSON linestring segments in OpenLayers?

texts:
I am rendering a vector linestring via GeoJSON on top of an OpenLayers map.  I'd like to alternate the styling between two slightly different shades of blue, such that the colors alternate from one line segment to the next.  Is there any way to do this, apart from using points instead of a linestring and constructing the line segments myself?

var path;

function drawPath(table, x, y, id) {
  var urlreq = "../grab.php?method=get_path&x=" + x + "&y=" + y + "&table=" + table +
                "&id=" + id;
  if( typeof(path) != 'undefined' ) { map.removeLayer(path); }
  path = new OpenLayers.Layer.Vector("Path", {
    projection: epsg900913,
    style: {
      strokeColor: "blue",
      strokeWidth: 3
    },
    strategies: [new OpenLayers.Strategy.Fixed()],
    protocol: new OpenLayers.Protocol.HTTP({
      url: urlreq,
      format: new OpenLayers.Format.GeoJSON()
    })
  });
  map.addLayer(path);
}


texts: Levels in NetCDF created in R with raster package are interpreted as time steps - create 4D rasterBrick

texts:
I can't seem to find a way to create an empty rasterBrick and then create a NetCDF file that has no time axis and only vertical levels. Here's my code

mylevels=10
s1 <- brick(nc=180, nr=90,nl=mylevels)
for (i in 1:mylevels) {

    s1[[i]]<-rasterize(... processing....
}

writeRaster(s1,'myfile,nc',format='CDF',varname='var1',varunit='kg')

Now when I open the file with any NetCDF software such as CDO or ferret my levels are read as time steps instead of vertical levels. I have tried to use zvar='var1' but it's not working. I am trying to avoid the NetCDF package because my data processing works just fine. I also can't find any information on how to create a rasterBrick with 4 dimensions. Is it possible?


texts: OGR Shapefile to CSV: Using Python to get WKT for multipart Polys

texts:
Using ogr2ogr and python, I'm running into a wall with trying to create a CSV with WKTs for a shapefile that contains multipart polygons. Currently this is this code I'm using (found it on internet):

#Purpose: To export a shapefile to WKT

from osgeo import ogr
import sys, os

input = ogr.Open(sys.argv[1])

layer_in = input.GetLayer()
layer_in.ResetReading()
feature_in = layer_in.GetNextFeature()

outfile = open(sys.argv[1] + ".wkt", "w")

while feature_in is not None:

    geom = feature_in.GetGeometryRef()
    geom_name = geom.GetGeometryName()

    outfile.write(str(geom)+ '\n')

    feature_in = layer_in.GetNextFeature()

So, do I need to use a different geometry container or do I need to evaluation each feature and if it's a multipart use a different set of geometry classes? I've looked around in the gdal/ogr documentation but I'm having a hard time reconciling it.


texts: "Geom-type or SRID not allowed" errors when importing sqlite layer to DB manager in QGIS?

texts: I am trying to import an sqlite file to the QGIS DB Manager with default import settings. While the layer displays fine in the graphical interface, the DB Manager only partially loads it, stopping after 2K of 500K rows, citing that the table "violates Geometry constraint [geom-type or SRID not allowed]". What am I doing wrong?


texts: GeoServer - labels not working correctly

texts:
Currently having a few issues with GeoServer not displaying our labels correctly. The code in question is here:

<Fill>
          <CssParameter name="fill">
            <ogc:Literal>#000000</ogc:Literal>
          </CssParameter>
          <CssParameter name="fill-opacity">
            <ogc:Literal>0.75</ogc:Literal>
          </CssParameter>
        </Fill>
        <VendorOption name="spaceAround">-1</VendorOption>
        <VendorOption name="conflictResolution">false</VendorOption>
        <VendorOption name="partials">true</VendorOption>
        <VendorOption name="goodnessOfFit">0</VendorOption>

      </TextSymbolizer>

    </Rule>
  </FeatureTypeStyle>

</UserStyle>

That's at the end of my style file. From what answers I've already looked at, this should be okay. However, labels appear and disappear at different zoom levels, as below:

I've turned guttering up to 100px, and also the metatiling to 1x1, but this is still happening. This does appear to have occurred a few other times on posts, but I think I've done everything they've done (as described above).


texts: QGIS Graphic modeler: Loading in all the layers

texts:
Currently I have special script/model, where I convert a geoJSON feature collection to some seperate shapefiles.
There was small pre step and prerequisite, that I want to eliminate, normally we load this GeoJSON file into QGIS by loading all the layers in by adding it as a datasource, either using the QGIS Browser or manage data sources (CTRL+L). The we specify the points, line and polygons that go trought the script.
You would think I could just select this file, as a vector layer and then use the "Load in Project Layer" algorithm to make it much more easier.
Sadly it will only load one of the features of the GeoJSON (for some unknown reasons this is only the points). I would use the "File/Map" Parameter but that sadly does not connect to my script at all because that parameter is not considered a layer file.
So how create an algorithm script do I do "Load in Project Layer" but with a feature collection, loading all the layers in the feature collection?
Edit: I have found a possible solution for now, however there is one limitation with it. I will post this solution later when I get the time.


texts: ST_Union fails with TopologyException despite valid polygons and using ST_SnapToGrid

texts:
I am working with Postgres 9.6 and PostGIS 2.2.2. These are my tables:

                                     Table "public.allparcels"
  Column  |            Type             |                        Modifiers
----------+-----------------------------+----------------------------------------------------------
 gid      | integer                     | not null default nextval('allparcels_gid_seq'::regclass)
 the_geog | geometry(MultiPolygon,4326) |
 county   | character varying           |

                                       Table "public.parcels"
    Column    |          Type           |                            Modifiers
--------------+-------------------------+------------------------------------------------------------------
 ogc_fid      | integer                 | not null default nextval('parcels_ogc_fid_seq'::regclass)
 wkb_geometry | geometry(Polygon,4326)  |
 county       | character varying       |

When I run this (via psycopg2):

INSERT INTO allparcels(county, the_geog)
SELECT 'Leeds', ST_MakeValid(ST_Union(wkb_geometry)) FROM parcels WHERE county='Leeds'

I get this error:

Traceback (most recent call last):
  File "combine.py", line 21, in <module>
    cursor.execute(q, (county, county))
psycopg2.InternalError: GEOSUnaryUnion: TopologyException: Input geom 1 is invalid: Self-intersection at or near point -1.394465365015132 53.741246412906918 at -1.394465365015132 53.741246412906918

I thought I must have invalid underlying geometries, but I can't see any:

SELECT * from parcels WHERE NOT ST_isvalid(wkb_geometry);

produces no results. As the query shows, I'm already running ST_MakeValid to make sure the unioned polygon is valid before inserting it.

What am I doing wrong?

UPDATE:

This is the latest query I have tried.  Debugging the query indicates that it's the ST_Union that fails. Decreasing the precision (even to 0.01) does not help:

INSERT INTO allparcels(county, the_geog)
SELECT 'Leeds', ST_MakeValid(ST_Union(ST_SnapToGrid(wkb_geometry, .000001)))
FROM parcels WHERE county='Leeds'


texts:
I need help in this question. To explain briefly -- the question is:

Whether a point with coordinates

lng:-80.087535
lat: 42.054246

lies over or under line with coordinates:

lng: -80.267831 ... -73.362579
lat: 42.050312

If it lays under (as Thomas said) then could you please show me where I've made a mistake in my coordinate system vision?

Update:
Sorry, my mind after hours of thinking is not working now. I haven't explained base of the question.
I've got two rectangles: the big one(R1) and a small one(R2). R2 is a part of a R1. Also I've got a point(P1) which is seems to be in R2, as MongoDB query says, and is missing in R1 when querying (as was explained in question - because of earth curvature).


texts: How do you use folders as input/output parameters

texts:
What I want to do is to go into a folder, find a layer file, pull it into an MXD, then export a jpeg in an output folder.  Everything is great except I don't know how to set the input/output folders.

state =  arcpy.GetParameterAsText(0)
In_Folder = arcpy.GetParameterAsText(1)
Out_Folder = arcpy.GetParameterAsText(2)

lyrfile_o = arcpy.mapping.Layer(In_Folder + str(state) + "_O.lyr")
lyrfile_d = arcpy.mapping.Layer(In_Folder + str(state) + "_D.lyr")
lyrfile_lane = arcpy.mapping.Layer(In_Folder + str(state) + "_Lanes_.lyr")
lyrfile_poly = arcpy.mapping.Layer(In_Folder + str(state) + "_Lanes_Polys1.lyr")

and I get an error like:

<type 'exceptions.ValueError'>: Object: CreateObject Layer invalid data source>

It's written in a text editor then brought into a toolbox and parameters set, both types as folders, the first one (state) is a string.

Any ideas?


texts: Adding Animation/Smooth Pan On Setting New Center on ArcGIS JS API Map

texts:
Can you please take a look at this demo and let me know how I can add smooth animation or Pan function to move to the destination?

As you can see I am trying to re-center the map center by

 map.centerAt(new Point(-78.542, 37.61, new SpatialReference({wkid: 4326})));

on a click event but this is jumping the map to the destination.

<!DOCTYPE html>

    <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
    <meta name="viewport" content="initial-scale=1, maximum-scale=1,user-scalable=no"/>
    <title>Simple Map</title>
    <link rel="stylesheet" href="https://js.arcgis.com/3.23/esri/css/esri.css">
    <link rel="stylesheet" href="https://maxcdn.bootstrapcdn.com/bootstrap/3.3.7/css/bootstrap.min.css" integrity="sha384-BVYiiSIFeK1dGmJRAkycuHAHRg32OmUcww7on3RYdg4Va+PmSTsz/K68vbdEjh4u" crossorigin="anonymous">

     <button type="button" class="btn btn-default">Go There</button>


texts: How to fix the "Snap tolerance" error message in QGis?

texts:
I am trying to use the "Delete part" tool on a shapefile. While I am starting to use it, an error message shows up : "Could not snap segment. Have you set the tolerance in Settings > Project Properties > General?". I tried to set different tolerances in the snapping options, the error still shows up and I can't delete any part of my polygon.
How to fix that?

I use QGis 1.8.0-Lisboa and Windows 7.

Thanks


texts: How to access ArcSDE metadata using OGR?

texts:
I'm trying to get a list of all layers in an ArcSDE database hosted in oracle

here is the command I'm using:

ogrinfo -al SDE:localhost,5151,sde,sde,sde

Note: I've replace the database, username and password values with sde

When I execute the above I get the following message:

FAILURE:
Unable to open datasource `SDE:localhost,5151,sde,sde,sde
Following drivers.
  -> OCI
  -> SOSI
  -> ESRI Shapefile
  -> MapInfo File
  -> UK .NTF
  -> SDTS
  -> TIGER
  -> S57
  -> DGN
  -> VRT
  -> REC
  -> Memory
  -> BNA
  -> CSV
  -> NAS
  -> GML
  -> GPX
  -> KML
  -> GeoJSON
  -> Interlis 1
  -> Interlis 2
  -> GMT
  -> GPKG
  -> SQLite
  -> ODBC
  -> WAsP
  -> PGeo
  -> MSSQLSpatial
  -> OGDI
  -> PostgreSQL
  -> MySQL
  -> PCIDSK
  -> OpenFileGDB
  -> XPlane
  -> AVCBin
  -> AVCE00
  -> DXF
  -> Geoconcept
  -> GeoRSS
  -> GPSTrackMaker
  -> VFK
  -> PGDump
  -> OSM
  -> GPSBabel
  -> SUA
  -> OpenAir
  -> PDS
  -> WFS
  -> HTF
  -> AeronavFAA
  -> Geomedia
  -> EDIGEO
  -> GFT
  -> GME
  -> SVG
  -> CouchDB
  -> Idrisi
  -> ARCGEN
  -> SEGUKOOA
  -> SEGY
  -> XLS
  -> ODS
  -> XLSX
  -> ElasticSearch
  -> PDF
  -> Walk
  -> CartoDB
  -> SXF

I'm assuming that since SDE isn't in the list that is the problem.

Any idea of how to add this support?

according to http://www.gdal.org/drv_sde.html

It isn't compiled by default and requires the ArcSDE drivers

I'm unsure how to install this correctly


texts: QGIS how to use the NNJoin Plugin algorithm within a Python script?

texts:
I have two layers: 1 point layer and 1 line layer. My goal is to calculate the nearest line to every point and I am doing this with the NNJoin Plugin. Now I want to do the same within a Python script, but I do not know how to access the NNJoin-Algorithm from here.

If this is not possible, are there any alternatives?

I have already read through the duplicate question but it doesn't quite answer my question.


texts: Why is my OnChangeFeature event not firing in my parcel fabric editor extension?

texts:
I have an edit extension that listens to events that happen during an editor session. I have wired it like so:

void Events_OnStartEditing()
        {

            //Since features of shapefiles, coverages, and so on, cannot be validated, ignore wiring
            //events for them.
            if (ArcMap.Editor.EditWorkspace.Type !=
                esriWorkspaceType.esriFileSystemWorkspace)
            {
                //Wire OnCreateFeature edit event.
                Events.OnCreateFeature += new IEditEvents_OnCreateFeatureEventHandler
                    (Events_OnCreateFeature);
                //Wire onChangeFeature edit event.
                Events.OnChangeFeature += new IEditEvents_OnChangeFeatureEventHandler
                    (Events_OnChangeFeature);

            }

        }

The OnStartEditing event fires just fine. I want to listen for when an attribute from the row changes, so I wired an OnChangeFeature event:

void Events_OnChangeFeature(ESRI.ArcGIS.Geodatabase.IObject obj)
        {
            MessageBox.Show("Event fired");
        }

Every time I go in and change a feature by editing the attribute table, this should fire. But it doesn't. I am manually editing features of the table (I'm not using a cursor to programatically update the features).

I will lastly mention that the feature class that I am editing is a subset of a parcel fabric. I'm not sure that this makes a difference as there is nothing in ESRI documentation that I have found that says that editing in a parcel fabric bypasses normal editing events, but I will mention it anyway.

EDIT: Ideas I have tried with no result:

Implementing IObjectClassInfo2 in my extension to make sure that no store methods are being bypassed:

public class ParcelEditorExtension : ESRI.ArcGIS.Desktop.AddIns.Extension, IObjectClassInfo2

        public Boolean CanBypassStoreMethod()
        {
            return false;
        }

        public Boolean CanBypassEditSession()
        {
            return false;
        }

Running the IWorkspaceEditControl.SetStoreEventsRequired() method to enforce this rule over the entire edit session.

IWorkspaceEditControl EditControl = ArcMap.Editor.EditWorkspace as IWorkspaceEditControl;
EditControl.SetStoreEventsRequired();

EDIT 2: It appears that the parcel fabric is, in fact, the issue. The event fires just fine when working on simple feature classes. Regardless, setting CanBypassStoreMethod() to false should fix this issue - is my implementation of this method correct?


texts: Using script to list .lyr files with broken links?

texts:
I would like to identify within a particular folder the .lyr files that have broken links. Some of these files are within multiple folders but all fall within one  folder. My programming knowledge is fairly novice so I may use confusing or incorrect terminology.

I have been following the Q&A below, and have got the script to successfully identify .mxd's in a test folder only containing one .mxd file containing a broken .lyr file. When trying the .lyr script given in the same thread on a test folder containing a broken .lyr file both IDLE and Arc Desktop (10.2) crash.

Use arcpy.mapping to list broken data layers?

successful broken .lyr in mxd script:

import arcpy, os
path = r"M_XXX"
for root, dirs, files in os.walk(path):
    for fileName in files:
        basename, extension = os.path.splitext(fileName)
        if extension == ".mxd":
            fullPath = os.path.join(root, fileName)
            mxd = arcpy.mapping.MapDocument(fullPath)
            brknList = arcpy.mapping.ListBrokenDataSources(mxd)
            for brknItem in brknList:
                print "MXD: " + fullPath
                print "\t" + brknItem.name

crash .lyr script:

import arcpy, os
path = r"M_XXX"
for root, dirs, files in os.walk(path):
    for fileName in files:
        basename, extension = os.path.splitext(fileName)
        if extension == ".lyr":
            fullPath = os.path.join(root, fileName)
            mxd = arcpy.mapping.MapDocument(fullPath)
            brknList = arcpy.mapping.ListBrokenDataSources(mxd)
            for brknItem in brknList:
                print "LYR: " + fullPath
                print "\t" + brknItem.name


texts: Access resolved/rendered dynamic text from ArcPy or ArcObjects?

texts:
Using ArcPy (or ArcObjects if essential), is there a way to retrieve the rendered text of a text element that includes Dynamic Text?

The text property of the TextElement includes the Dynamic Text tags, but I would like to access the resolved text that will be displayed on the map.

For example for a title page at the moment:

myTextElement.text = "Country Map - Name: <dyn type="page" property="name"/>
                     - Code: <dyn type="page" property="attribute" field="CNTR_ID"/>"

But I'd like to get the resolved dynamic text, in this example from a data driven page:

Country Map - Name: Zimbabwe - Code: ZW

This is to be able to export metadata for maps produced using data driven pages.  At the moment without a method to access the resolved value of the dynamic text the only solutions we can think of are to either;

Write a class to parse and resolve the dynamic text.

This wouldn't be too difficult to do as the datadriven page class includes a property to the index layer, but I expect running an attribute cursor for each page would slow exports down.

Save the MXD for ArcMap 9.3 which doesn't support dynamic text (this strips the dynamic text from the layout and replaces it with it's values).

While this would work it seems a long way around, as after saving a copy the map document would then need to be opened to retrieve the text element value.


texts: Get convex hulls for all rasters in raster_columns

texts:
how to search the polygons of the images that I have in the Database?
I want execute this function ST_AsText(ST_ConvexHull(rast)) for every "r_table_name"

SELECT
    r_table_name
FROM
    raster_columns

SELECT
   ST_AsText(ST_ConvexHull(rast))
FROM
    NDVI_3_fev_18


texts: Extract specific values from raster in QGIS

texts:
in QGIS 2.18.3, Windows 10, I've got a raster file in which each color defines a type of crop. See image:

Each cell value represents a type of crop. What I need to do is to extract to a new file, only the cells containing a certain value (e.g., the value 24 which corresponds to the brown cells) Ultimately, the objective is to calculate the surface of each crop within a given area.

Is there any process that can help me with this analysis?


texts: Converting a .dxf to a .shp containing a field with height information in QGIS?

texts:
I received two .dxf for my project that represent my area of interest: one of them contains buildings, the other one terrain elevation (they both are 3D). I can successfully convert them to shapefiles using the Dxf2Shp plugin. They also display nicely compare to sample below:

However, I loose the height information if I do so. The only field I get in the attribute table is "myid", that simply counts the polygons.
Is there a way to convert the files to .shp creating a height field in the attribute tables, that takes the z values from the .dxf? This would already be fine, as the projection is metric.


texts: Matching two point layers (with different gridding) to common grid?

texts:
I have two different point layers, where points represent centroids of areas covering a stream segment. One layer (lets call it layer A) is coarser, with fewer point (each point represents a larger area), then layer B. There is no commonality between how the points are distributed, except for the boundaries. Each layer contains multiple attributes (ex. depth, velocity, areas). I want to create a common grid, but I'm hesitant about rasterizing this data, as I need to retain the multiple attributes.

I have tried converting each set of points to voronoi diagrams, then spatially joining layer A to layer B, and converting back to points, but there have been problems retaining area values required.

Also I have joined layer B's voronoi areas to layer A's points. This gives me layer B's areas in terms of layer A's points, but leaves holes in the coverage.

Ideally, I'd like layer A to fit within layer B, or have a new mesh fitted to the stream segment that both these layers can be applied to.

Does anyone have suggestions for optimizing this process, or ideas I've failed to see/am not familiar with?

I've included pictures of the point layers, as well as a visual goal.

Layer A picture

Layer B picture

Solution?

I'm using ArcGIS 10.3 with Spatial Analyst, and QGIS 2.4 Chugiak.


texts: Infrastructure networks in GIS modelling

texts:
I am currently working on a GIS-based project which aims to model the interdependency of existing infrastructure networks, e.g. transport, oil, pipeline networks, etc. The network interdependency include the interactions between networks in terms of physical or  geographical respective. This is quit important as it might lead to cascading effect among these networks.

We would like to develop a modelling platform/software with GIS function so the advantages of geographical relationships can be fully explored. We have already started some works based on a few open-source libraries and software. However, I hope I am open-minded, and am wondering if there are any existing software tools/publication related to this aspect?


texts: Cropping vector map into tiles with QGIS - how to get around shapefiles which have stray points on boundary of tiles?

texts:
I am trying to use QGIS to generate a set of tiles in Shapefile format from a vector coastline map.

My script is correctly generating the 5x5 degree tiles I want, however a few tiles clearly have some data missing. Investigating these tiles we noticed that they have lone points on the boundary of the tile (e.g. from a peak in the neighbouring tile which just clips into the tile by one point). This is considered a failure and the entire feature (which connects to that point) is removed from the tile.

Output from the command

"gdalogr:clipvectorsbyextent', layer, extent, "-nlt " + ltype + " -skipfailures -overwrite", file"

shows we are getting this type of error:

Attempt to write non-linestring (GEOMETRYCOLLECTION) geometry to ARC type shapefile

This is due to the presence of Linestrings and a stray Point.

I've attached some images which demonstrate this. The lighter green line is the entire map we are cropping from, and the dark green shows the lines that have been kept when generating that tile. If you zoom in on the line just touching the left side of the tile you can see there is one point of that feature overlapping into our tile.

What would be the best way to get the map tiles in shapefile format and not have any missing data? I also need to have only one shapefile per tile.


texts: Align in print composer not working in QGIS?

texts:
When i try to use the Align tool in Print composer to align either labels, arrows or maps, nothing happens. I've tried restarting QGIS, but it seems to be the same issue in every project. What am I doing wrong here? Or could this be a bug?

Using QGIS 2.10.1 in Windows 7.


texts: GDAL Virtual File Systems Python Opening files in Network Drive Issue

texts:
I'm trying to open a .tar file using /vsitar/
The file to be opened is located at: //my_network_drive/scratch/my_file.tar
In Python I have the following snippet:
driver = ogr.GetDriverByName('ESRI Shapefile')
driver.Open(/vsitar//my_network_drive/scratch/my_file.tar/my_file.shp`

This does not seem to work even though the file exists.
However, if I were to map the network drive say to Z
driver = ogr.GetDriverByName('ESRI Shapefile')
driver.Open(/vsitar/Z:/scratch/my_file.tar/my_file.shp`)

This does work.
The documentation seems to show that what I want to achieve is possible
/vsitar//home/even/my.tar/subdir/my.tif # (absolute path to the .tar)
https://gdal.org/user/virtual_file_systems.html#vsis3
Does anyone know a way to fix this without having to specifically map a drive?


texts:
I have the following model which uses two submodels that have iterators inside them. I need to spatially join the outputs of these models. Currently, it will completely run the First_Submodel, overwriting all the outputs in the process. Then it will completely run the Second_Submodel. I need it to run the First_Submodel and create Zone1. Then the Second_Submodel will start and create Zone2. Then it will spatially join the two outputs and that will be the first result. Then it will again start the First_Submodel, create Zone1(2), start the Second_Submodel, create Zone2(2), spatially join them, and create the second result. And so on.
Somewhere in the model, there should probably be inline variable substations and preconditions, but I am not sure where. Is it even possible to make the ModelBuilder run that way?


texts: I am currently using the QGIS C++ libraries to display a map. I am new to QGIS, and wantd to know what the difference between map and layer coordinates is? From what I understand, map gives the coordinates as per the CRS of the canvas. Is that correct?


texts: photo2shape loading error classFactory()

texts:
Python error: Couldn't load plugin photo2shape due to an error when calling its classFactory() method See message log (Python Error) for more details.

Just updated to the newest version of QGIS and upon trying to use photo2shape when I got this message. Does anyone have an idea s to how to solve the issue?

I'm using version 2.16.2


texts:
I created a custom web TileServer app that serves map tiles using a specific url in the form

localhost/2/{z}/{x}/{y}.png

z, x, y are tile coordinates (like for instance in OpenLayers or Google)

How can I add a layer from this tileserver to QGIS?


texts: I have aerial photos that are .jp2 format and are compiled into a .dbf. In ArcGIS I can simple add the .dbf file and all the aerials automatically appear. Can this be done in QGIS, rather than insert each individual .jp2 file?


texts:
I'm using the approach described here to change LiDAR point cloud projection in R:
Can I re-project an LAS file in LidR
But, what about lidR catalog reprojection?
If I used the same reproducible dataset than in the post Can I re-project an LAS file in LidR:
library(lidR)
LASfile <- system.file("extdata", "Megaplot.laz", package="lidR")

# las approach
las = readLAS(LASfile)
las
las2 = spTransform(las, sp::CRS(SRS_string = "EPSG:26918"))
class(las2)
# [1] "LAS"
# attr(,"package")
# [1] "lidR"

# lascatalog approach
ctg <- catalog(LASfile)
ctg2 <- spTransform(ctg, sp::CRS(SRS_string = "EPSG:26918"))
class(ctg2)
# [1] "SpatialPolygonsDataFrame"
# attr(,"package")
# [1] "sp"

The result when using the same command with a lascatalog is a SpatialPolygonsDataFrame and not a lascatalog... Is there a way to reproject lascatalog while remaining in R?
I don't like the las2las solution because you may only work with projection already in pre encoded by Martin within the tool (or maybe I'm missing something with the use or las2las)


texts: Calculating a specific distance from points along a network - GRASS

texts:
I'm currently using the v.net.iso function in GRASS (no access to network analyst in ArcMap). Does anybody know of a way to export the result to shapefile so that you have an individual line feature(s) for each point?
Currently, one can only see what areas along the network that are within that specific distance from the nearest point as opposed to being able to select an actual line segment of specified distance from each point.

The current result does not show you an actual specified distance path, but rather where along the network a segment of the network falls within that distance (difficult when there are overlaps - no way of seeing where these distances actually fall).


texts: What to use as split field or target workspace with Split tool?

texts:
I'm trying to use the split tool in arcgis 10.1 but I dont know what I should use as the split field or my target workspace.

W:\S&P\s&p techs\Emily\Geodatabase\Geodatabase.gdb
this is what I have as my target workspace but it says it's empty.


texts: Formatting SQL expression for calculate field (management) tool?

texts:
I am trying to take the population census tract and calculate the percentage of population by each census tract in the field (Populati_2). To do this I'm creating a new field with the total sum of the population (pop).

To get the percentage I'm usingthe following expression on arcgis [Populati_2] / [Pop] *100. It works fine in arcgis but when I'm using it to calculate field management or python it just doesn't work.

I'm not sure how to format the expression, I have tried every single combination and it's just not working.

Here's part of the script:

arcpy.AddField_management(r"Results\popjoin.dbf", "Pop", "LONG", "", "", "5")
arcpy.AddField_management(r"Results\popjoin.dbf", "Pop2", "DOUBLE", "", "", "5")
arcpy.CalculateField_management(r"Results\popjoin.dbf", "Pop", summed_total)
arcpy.CalculateField_management(r"Results\popjoin.dbf", "Pop2", '''Populati_2] / [Pop] *100''')

This one gives me errors, but the following does not give me any errors but I get a 0 in every row instead of the calculation.

expression = "!Populati_2! / !Pop! * 100"

arcpy.CalculateField_management(r"Results\popjoin.dbf", "Pop2", expression, "PYTHON_9.3")


texts: ComException -2147467259 Arcgis 10 VB.net +Geoprocessing

texts:
I have a problem when i execute GeoProcessing in ArcEngine 10 with VB.net

i had create the toolbox : Models.tbx
and  Model Process : disolveRoute with 2 parameters

My Function is :

 Sub diso()
    Dim dissolveTable As ITable = CType(featureLayerRoute, ITable)
    ' Create geoprocessor. Overwrite true will replace existing output.
    Dim gp As IGeoProcessor2 = New ESRI.ArcGIS.Geoprocessing.GeoProcessor()
    Dim trac As ITrackCancel = New CancelTracker
    Dim Result As ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult2
    Dim parameters As IVariantArray = New VarArrayClass()
    gp.SetEnvironmentValue("workspace", "Database Connections\IFE connexion.sde")
    gp.OverwriteOutput = True
    gp.AddToolbox("data\Models.tbx")
    parameters.Add(dissolveTable)
    parameters.Add("typeroute")
    Try
        Result = CType(gp.Execute("disolveRoute", parameters, trac), ESRI.ArcGIS.Geoprocessing.IGeoProcessorResult2)
    Catch ex As COMException
        Console.WriteLine(ex.ErrorCode)
    End Try

End Sub

For the inpute feature for disolve --> featureLayerRoute :

  Dim Sp As IWorkspace = ConnectToTransactionalVersion("localhost", "sde:postgresql:localhost", "sde", "....", "....")
        featureLayerRoute.FeatureClass = featureWorkspace.OpenFeatureClass("piste_route")
        Dim layerrouteAs ILayer = featureLayerDroit
        layerroute.Name = featureLayerRoute.FeatureClass.AliasName

but when i execute this function it returned a Code ERROR -2147467259.. Please it's one Week i'm looking to solve it but Nothing :'(

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


texts: What are the units in the network analysis

texts:
.

I'm working on a project that is checking to see if distance is related to colonization by invasive species based on human activity.  In order to do that, I must calculate the distances between colonized lakes and uncolonized lakes in a variety of forms.  I have done this using QGIS 3.12 by using the built in network analysis run point to point, and then the QNEAT3 plugin using the OD Matrix from Layers as Table (m:n).  Both result in very similar values, but these values don't have units and I am unable to figure out what they are.  The layers have units of degrees, but can't be degrees.  Google and StackExchange searching has not lead me to any answers, and neither did reading the manual for QGIS or QNEAT3.  I could guess at meters,but a distance of 190-191 km is way to high. If it is in feet,its way to low. It's actually only about 11 miles or so between the points.  I could just go with the "cost" values, but I need actual units, not values.

Doing one or two of these wouldn't be bad, but I have potentially thousands that will need to be calculated set over many states, and I don't have 5 years to do that on google.

If anyone could let me know either what the units are, or how to convert them to km or miles, or both, please let me know.

Thanks,

Greg W.


texts:
I am working on a twelfth-century description of South Italy (Calabria, Apulia and Campania) which describes rivers that are now dried up.
I am attempting to create maps which reflect the description from this twelfth-century text in ArcMap.
When I look in Google Earth I can see these dried up rivers. Is there a place to find data on all the rivers in Italy, including dried up rivers?


texts: Adding labels to contours created with the rendering option of rasters in QGIS 3.16

texts: Is it possible to add labels to raster contours created using this option in the "Render Type" under "Raster symbology"?


texts: Closing gaps between two rasters in QGIS

texts:
I have two rasters which are supposed to be perfectly aligned but they are not. There are 'gaps' between them as you can see in the picture (the white stripes are gaps):

How can I align them so that they fit?
My first thought was to use a resample tool or Align Raster but it didn't bring the right result.


texts: How to convert meters in height to a z axis/elevation

texts: I have an dataset with buildings and their height in meters. The epsg is 28992. is there a way to convert their height in meters to a z axis that represents the heights of the buildings in meters in python.


texts:
I am currently using the OpenLayers plugin to import Google Maps as the base layer for a heat map. The map is of population by census tract. The more populous census tracts are red and less populous green.

This is working well but the OpenLayers is completely flat and 0% transparent. This means I have to have my coloured census tracts over the top which obscures the features of the map.

In products like MapPoint you are able to have the background of the map at the very back, then the coloured data and then the labels (cities, highways etc.) at the top.

Does anyone know how I can do that?

Another option would be a transparent map layer which I could then put over the top of my coloured census tracts.

Any help would be much appreciated.


texts: How to do effective spatial joins using spark v2.0 and python?

texts:
I need to do spatial joins and KNN joins on big geolocalised dataset.

They are stored as pySpark RDDs.

I've been looking for libraries to do so, but couldn't find any that fits my needs:

compatible with Spark 2
compatible with pySpark
Allowing to do fast spatial joins

Have anybody succeed to do geo-analysis with pySpark ?

Would using Spark in combination with shapely functions or geoPandas dataframes be a decent solution ?

PS: Libraries I looked :

Magellan => not compatible with Spark2
SpatialSpark => neither
GeoSpark => not compatible with pySpark
LocationSpark => neither


texts: Check if a single feature is completely within another (arcpy)

texts:
I have a feature class that contains one polygon feature, and I need to check if said feature lies completely within another polygon feature.

Used selectLayerByLocation and then retrieved the count (arcpy.GetCount_management) for some other intersection based checks that need to be done; in this case however there is only ever one feature in the feature class and retrieving the count will only return one (regardless of if the polygon is fully contained or not).

Any ideas?


texts: how to add timeline controls to cartodb torque heatmap

texts: I figured out how to add multiple torque heatmap layers to a map using  cartodb.js but I can't figure out how to add a timeline control. I know it's probably something than can be done with one line but I can't find it in the documentation and resources.


texts: Can ArcREST be used with a local ArcGIS Server?

texts:
I'm looking at using ArcREST to insert rows in a featureservice on a local ArcGIS Server.

All the sample that I've seen seem to be connecting to a ArcGIS Online.

If I look at the Documentation for the Security Handler Helper module,
I can't see where to specify the url of the local server.

Is there a Sample for inserting rows in the featureservice of a local ArcGIS Server?


texts: Getting broken up area of smaller polygon that overlaps two larger polygons using ArcMap

texts:
I am using ArcMap.
I have two polygons. One is precinct which there are around 45 and i have census tracts which there are around 1,0000.
I want to create total population within each precinct. Though many census tract can be easily added since they fit in, not all are completely contained within one precinct. Some census tracts overlap two or even three precincts.
If they fit in completely the population value can be added no problem, however i want to split the population by percent of the census tract within a precinct.
For example if i have one census tract at a value of 100 and it crosses precinct 1 and 2, at 25% in precinct 1 and 75% in precinct 2 then I want a value of 25 allocated to precinct 1 and 75 in precinct 2.


texts: Is there an itemSelectionChanged Signal in Attribute Table?

texts:
The QgsVectorLayer class emits a selectionChanged() signal . If a QGIS user for example selects a record in the Attribute Table of a Vector Layer - a selectionChanged() signal is emitted an can be connected to an appropriate Slot inside the python code of a plugin. This signal is emitted every time the selection changes. I would need to realize a similar approach, but from a different direction. A signal should only be emitted when a user interaction occurs.

For example:

The user opens an attribute Table (of a vector layer) and selects one ore more records in this Table View with a mouse click. Is it possible to receive this click or select signal within a python plugin?

Does the Attribute TableView  in QGIS itself emits a kind of itemSelectionChanged Signal whenever the selection of  Attributes in the Table View changes? Which classes could be used to get this signal? I couldn't find any Information concerning my problem in the QGIS-API Documentation.


texts:
Using Putting content of popup into modal using OpenLayers? I am now able to put the content of a popup into a  modal of Bootstrap 4.

...

var overlay = new ol.Overlay({
  element:  document.getElementById('myModal'),
  autoPan: true,
  autoPanAnimation: {
    duration: 250
  }
});

var map = new ol.Map({
  controls: ol.control.defaults().extend([
    new ol.control.ScaleLine({className: 'ol-scale-line', target: document.getElementById('scale-line')}),
  ]),
  layers: [
    OpenStreetMap,
    istat,
  ],
  overlays: [overlay],
  target: 'map',
    render: 'canvas',
  view: new ol.View({
    extent: boundary,
    center: center,
    zoom: 2,
        maxZoom: 20,
        minZoom: 14
  })
});

map.on('singleclick', function(evt) {
  clickFeature = map.forEachFeatureAtPixel(evt.pixel, function(feature, layer) {
    return feature;
  });
  if (clickFeature) {
    var element = overlay.getElement();
    overlay.setPosition(evt.coordinate);
    $(element).modal('show');
  }
});

  @TomazicM provided part of this code

The problem now is the animation of the map that is activated after the click on the geometry. If I click on a geometry the modal opens and the map moving down. In other words happens a thing similar of the use the pan function for navigate the map.

How I can remove this?


texts: Printing CJK characters in title using ArcGIS API for JavaScript?

texts:
I am trying to print Chinese characters title using ArcGIS javascript (3.10); however, the characters becomes a square block. I am able to see the characters in the JSON passed to the print task but the output is incorrect.

printTemplate.layoutOptions = {
  titleText = "some chinese characters"
}

Is there a way to set the print task to use Unicode characters?


texts:
I have a street map, from Durango, Mexico, and I want to print it. But QGIS show the next error "EXCEPTION bad allocation". I'm going to print it in 1.50 x 1.06 mts sheet.

What can I do?


texts:
I have a spatial database in Sql server and I am working in asp.net mvc. My project must show stored shape in web page that uses openlayers js and user must see and edit shapes and send the result to server. I want to know how can I send data to my map and receive new data from it and save data to database.

Is there any library that work for me? Is is not important that it must be open source

thanks a lot


texts:
Does anyone know if Australian topo map GeoTIFFs are made available anywhere (preferably for free)?

LINZ makes them available in New Zealand - http://www.linz.govt.nz/topography/topo-maps/index.aspx, and I was hoping a similar service was available for Australia.