texts: Custom Tools not Saving in Toolbar on Startup/ Normal.mxt?

texts:
I am trying to save some tools that I have created into the default toolbar layout that is present when ArcMap 10.1 is launched.

It is possible to change the arrangement of toolbars, add new ones and even add in existing tools, but any custom tools I set in the commands category "geoprocessing tools" are gone and do not appear in the toolbar. How can I save these custom tools to the default layout that appears on startup?

Update:

Saving a document with the custom tools added and reopening it preserves the tools, but opening a template (10.1 mxd) containing the tools, then opening and older document will cause the tools to disappear. Oddly enough, any other changes made after the custom tools were added are preseved.

Update (03/18/2014):

I have upgraded to 10.2 and this has still not been fixed. I did not see this issue in the 10.21 changelog so I doubt it will address it. For those interested in this issue, I have started a forum post on the ESRI forums here.


texts:
I need to reclassify unsigned integer raster data.

I ran Calculate Statistics on the Raster, tried Reclassify, but keep getting  Error 999999: Error Executing Function. The description for the error is "This is a generic error for which the cause does not have a specific error ID."

I have 10 categories and need 2-10 reclassified to "0." Category 1 remains "1"

Here's my environments (I'm using ARCMAP 10.3.1)

What might the problem might be?


texts: how to use union tool for merging more than two layers?

texts: I want to combine different layers into one but I still want to identify different layers in the new combined layer. Union is best for this, but it doesn't allow me to combine more than two layers. I am using arcgis 10.1, with site license.


texts:
I'm trying to use Bruce Harold's Concave Hull Estimator Tool and I'm getting some errors.

When trying to use the tool on ArcMap, I get this:

There seem to be no object that I can use the tool on, even though what I have loaded is basically points with two attributes:

Data Type:  Shapefile Feature Class
Shapefile:  pts.shp
Geometry Type:  Point
Coordinates have Z values:  No
Coordinates have measures:  No

Projected Coordinate System:    WGS_1984_Web_Mercator_Auxiliary_Sphere
Projection: Mercator_Auxiliary_Sphere
False_Easting:  0.00000000
False_Northing: 0.00000000
Central_Meridian:   0.00000000
Standard_Parallel_1:    0.00000000
Auxiliary_Sphere_Type:  0.00000000
Linear Unit:    Meter

Geographic Coordinate System:   GCS_WGS_1984
Datum:  D_WGS_1984
Prime Meridian:     Greenwich
Angular Unit:   Degree

Is there some other specific change I need to make to use this tool? Is there a better tool for concave hull that's easy to use??


texts: Lidar and infrared data to detect petroleum releases

texts:
Would this combo detect differences in soil temperature in proximity to the release due to the higher temperature of pipeline fluids or would the surface temperatures be affected by the volatilization of the released fluids?

I know this technique is being used currently and I am just curious about the process by which it would be able to determine a release point.


texts: R rasterEngine Parallel Processing, serialized and unserialized memory errors

texts:
I am trying to perform some raster math and wanted to try parallel processing in R using rasterEngine.
I am able to run the script for a few minutes and watch my CPU usage increase substantially, but soon an error occurs and processing stops.
Here is my code
library(raster)
library(spatial.tools)
library(doParallel)

cpus = 24 # The number of nodes/cores to use in the cluster
cl <- makeCluster(spec = cpus, type = 'PSOCK', methods = FALSE, outfile="E:/Heatload_10m/cluster.txt")
registerDoParallel(cl)

l = raster('E:/Heatload_10m/Lat_Radian.img')
s = raster('E:/Heatload_10m/R1_10m_WGS84_Slope.img')
f= raster('E:/Heatload_10m/Folded_Aspect.img')

heatload <- function(lat,slope,folded){
    hl= exp(-1.467 + 1.582*cos(lat)*cos(slope)-1*cos(folded)*sin(lat)*sin(folded)-0.262*sin(lat)*sin(slope)+0.607*sin(folded)*sin(slope))
    return (hl)
}

print('Starting Heatload Calc')
final_heatload <- rasterEngine(lat=l,slope=s,folded=f, fun=heatload, overwrite=TRUE, processing_unit = "chunk", chunk_format = "raster" setMinMax =TRUE, )
stopCluster(cl)
registerDoSEQ()

Here are the error messages:

The instruction at 0x61b83401 referenced memory at 0xd120c1d0. The memory could not be written.

Error in unserialize(socklist[[n]]) : error reading from connection
Error in serialize(data, node$con) : error writing to connection


texts: Arcobjects not drawing GraphicContainer if layers are in the TOC

texts:
Im trying to draw some graphic labels in Layout View for every row in a feature class chosen by a user in a combobox. I am just in testing out some code to draw the label graphic and it works fine, but as soon as I add layers to my TOC and try to click a button to generate the graphic it does not work and does not draw. Is there something I am missing?

        IMxDocument pMxDoc = ArcMap.Application.Document as IMxDocument;
        IActiveView activeView = pMxDoc.ActiveView as IActiveView;
        ITextElement pTextElement;
        IElement pElement;
        ICallout pCallout;
        IFormattedTextSymbol pTextSymbol;
        IGraphicsContainer pGraphicsContainer;
        IPoint anchorPoint = new PointClass();
        anchorPoint.PutCoords(50, 50);
        pTextElement = new TextElementClass();
        pElement = pTextElement as IElement;
        pTextElement.Text = "<UND><BOL>As-Cleared</BOL></UND>" + "\r\n" + " 3.20ha" + "\r\n" + "Partial Cleared";

        ESRI.ArcGIS.Display.IScreenDisplay screenDisplay = activeView.ScreenDisplay;
        ESRI.ArcGIS.Display.IDisplayTransformation displayTransformation = screenDisplay.DisplayTransformation;

        int x, y;
        displayTransformation.FromMapPoint(anchorPoint, out x, out y);
        IPoint messagePoint = displayTransformation.ToMapPoint(x + 120, y - 180);
        pElement.Geometry = messagePoint;

        pTextSymbol = new ESRI.ArcGIS.Display.TextSymbol();
        pCallout = new BalloonCallout();
        pTextSymbol.Background = pCallout as ITextBackground;
        pCallout.AnchorPoint = anchorPoint;
        pTextElement.Symbol = pTextSymbol;

        pGraphicsContainer = activeView as IGraphicsContainer;
        pGraphicsContainer.AddElement(pElement, 0);
        pElement.Activate(activeView.ScreenDisplay);
        activeView.PartialRefresh(esriViewDrawPhase.esriViewGraphics, pElement, null);


texts: Visualizing shapefiles on web using SHINY package of R with PostGIS?

texts:
I am using PostGIS database to store my shapefiles. Now I want to read those files through shiny package of R.

How can I visualize those files on web using shiny package of R and perform statistical analysis on those shape files?


texts: Effectively applying cpt-city colour ramps to elevation in QGIS3

texts:
There are lots of nice topography and bathymetry colour ramps provided by cpt-city in QGIS. However, I cannot find a straightforward way of specifying the transition from blue (below sea level) to greens/browns (above sea level) at (or in the region) of 0. Instead, I get something like this:
I could manually edit the values but this seems cumbersome especially over a large range.  Currently, the best option appears to be the one proposed here (split the data into 2, 1 above sea level, one below sea level):
Bathymetric-topographic DEM color relief?

Is there a way to set a certain value, e.g. 0, and then force the recalculation of the others? At the moment if I change the blue/green value to 0, it just reorders the colours. Am I missing something obvious in how these colour ramps should be applied?


texts: Mapbox Studio: How to upload glyphs or sprite for custom markers

texts: I am wondering how I can view the available or upload new glyphs or sprites into Mapbox Studio. I want to style some markers with custom icons and can't figure out where the default glyphs are coming from and how i can change that in Mapbox Studio.


texts: Using Python comtypes and ArcObjects to change dataframe background color in ArcMap?

texts:
I'm simply trying to programmatically change the background colour of a Dataframe in ArcMap.

I followed the beginning of this solution to access ArcObjects in Python:

import arcpy
import os

from snippets102 import *
from comtypes.client import GetModule, CreateObject

import comtypes.gen.esriFramework as esriFramework
import comtypes.gen.esriArcMapUI as esriArcMapUI
import comtypes.gen.esriCarto as esriCarto
import comtypes.gen.esriDisplay as esriDisplay

mapFile = "C:/MapMaps/myMap.mxd"

# Get map document object
pMapDoc = CreateObject(esriCarto.MapDocument, interface=esriCarto.IMapDocument)
pMapDoc.Open(mapFile)

From the Map Document object, I then tried to get the Active View, the Graphics Container (I don't entirely understand what that is), the Map, the Map Frame, the Background, the SymbolBackground, the Fill Symbol and finally the Color.

# Get Active View
pageLayoutActiveView = CType(pMapDoc.PageLayout, esriCarto.IActiveView)

# Get Graphics Container
graphicsContainer = CType(pageLayoutActiveView, esriCarto.IGraphicsContainer)

# Get Map
map = pMapDoc.Map[0]

# Get Map Frame
mapFrame = graphicsContainer.FindFrame(map)

# Get Background and Symbol Background
bgSymb = CType(mapFrame.Background, esriCarto.ISymbolBackground)

# Finally get to Fill Symbol and Color
print bgSymb.FillSymbol.Color.RGB

However, from here it doesn't seem I'm able to actually modify any of the Background symbol's properties

# Current color settings
print bgSymb.FillSymbol.Color.RGB
print bgSymb.FillSymbol.Color.GetCIELAB()

# Try to change...
bgSymb.FillSymbol.Color.RGB = 12345678
bgSymb.FillSymbol.Color.SetCIELAB(50, 10, -10)

# But settings are same as before!
print bgSymb.FillSymbol.Color.RGB
print bgSymb.FillSymbol.Color.GetCIELAB()

Any suggestions? I also tried creating an RGBColor object and doing bgSymb.FillSymbol.Color = newRGBcol but that did not work either.


texts: QGIS creating shapefile of several points out of two Excel sheets

texts:
I am about to create a QGIS file on which it should be visible which companies are working together. As source I have got two Excel tables.

Table = Distribution Centers
Table = Point of sales

I like to connect The point of sales (always several) with the associated distribution center. In the tables of the distribution centers you have the ID, and some (x/y) coordinates. In the Point of sale table you can find a ID, also (x/y) coordinates + a ID and which allows you to connect the Point of sales and the distribution centers (the second ID is the same like any of the distribution centers).

I want to separate the group of the distribution center and the related point of sales as one shapefile without searching manually in my Excel tables and creating 1 shape manually for every group...

I hope you will be able to understand my problem although my English is horrible.

My Lines are like this...

TDLinx-DC|TDLinx_-POS|Geometrie
00421362|00421362|LINESTRING (2.30567 48.7321)
00421362|00611434|LINESTRING (0.89378 46.4288)
00421362|00409507|LINESTRING (6.45527 46.359952)

Now it is like this. Might it be impossible to have a mixed point AND polylines layer?

TDLinx-DC|TDLinx-POS|Geo
00421102|00617429|LINESTRING (6.06043 48.70638, 7.34101 47.7495)
00421102|00256888|LINESTRING (6.06043 48.70638, 6.1812365 48.6876425)
00421102|00264459|LINESTRING (6.06043 48.70638, 5.0379298 47.3211782)
00421102|00264472|LINESTRING (6.06043 48.70638, 6.4946471 48.591343)
00421102|00264483|LINESTRING (6.06043 48.70638, 6.1618809 48.6763621)
00421102|00264501|LINESTRING (6.06043 48.70638, 7.7577037 48.5677807)
00421102|00264502|LINESTRING (6.06043 48.70638, 5.5904536 47.4484859)
00421102|00264503|LINESTRING (6.06043 48.70638, 6.1550176 47.6202079)
00421102|00264758|LINESTRING (6.06043 48.70638, 6.4474578 48.1758391)
00421102|00269360|LINESTRING (6.06043 48.70638, 6.0215766 47.2392361)
00421102|00269382|LINESTRING (6.06043 48.70638, 5.384374 49.1602425)
00421102|00269399|LINESTRING (6.06043 48.70638, 7.3564671 48.0791666)
00421102|00269478|LINESTRING (6.06043 48.70638, 6.8588071 47.6378244)
00421102|00429683|LINESTRING (6.06043 48.70638, 7.7465095 48.5839322)
00421102|00421102|MULTIPOINT (6.06043 48.70638, 7.34101 47.7495, 6.1812365 48.6876425)


texts: Extract ordered point coordinates from polygon layer via Virtual Layer or PostGIS View

texts:
I have polygon layer and I need to create table with XY coordinates for each polygon staring from farthest to west (lowest value of X).
Now I have Virtual Layer with

SELECT ST_DissolvePoints(geometry) AS geom, id
FROM "poly"

next I use processing to change multipart to singlepart layer I export them to CSV with XY coordinates and then I order them manually in excel.

Preferred solution is to have virtual layer with already ordered points but if solution with POSTGIS view will be good enough.


texts: Is it possible to do heat lines on MapServer similar to the following image. The datasource will be PostGIS.


texts: Orientation in 3D coordinate space of a pixel in georeferenced video imagery

texts:
I am trying to find (a) method(s) to extract the position in 3D coordinate space of a bounding box (detected object in 2D video imagery) having knowledge the below:

I have multiple images with slightly different angles and distances of the bounding box (detected object).
The multiple images are video frames that are products of 3 cameras.
I know the relative position between the cameras and a GPS receiver.
Every frame - image is "georeferenced" having knowledge about the exact position of the GPS.

In literature I am searching for stereosocpic and multiview techniques but is seems that I cannot find a suitable one.
Any help will be really useful..any idea "in a nutshell" will appreciated.

Thanks!


texts: Rasterize feature class in geodatabase using gdal.Rasterize (Python API)

texts:
I am trying to rasterize a polygon data set (feature class within a file geodatabase) using the GDAL Python API. I am using gdal.Rasterize (Actually, it would be very interesting to know what the difference between gdal.Rasterize and gdal.RasterizeLayer is. It doesn't seem to be explained anywhere).
This is the piece of code I am trying with:
from osgeo import gdal
from osgeo import ogr

#Geodatabase and Polygon Feature Class name (polygons to be rasterized)
path_gdb = r'.......\Flood_Defences.gdb'
fc = r'FloodDefences'

#Reading polygons with ogr
driver1 = ogr.GetDriverByName("OpenFileGDB")
gdb = driver1.Open(path_gdb, 0)
flood_def = gdb.GetLayer(fc)
sr = flood_def.GetSpatialRef()

#Extracting information from template raster
template_path = r'.........\template.tif'
raster = gdal.Open(template_path)
gt = raster.GetGeoTransform()
ndv = raster.GetRasterBand(1).GetNoDataValue()
proj = raster.GetProjection()
cols = raster.RasterXSize
rows = raster.RasterYSize
count = raster.RasterCount
dtype = gdal.GetDataTypeName(raster.GetRasterBand(1).DataType)

#Create object to store output
out_path = r'..........\rasterized_polygons.tif'
driver2 = gdal.GetDriverByName('GTiff')
output = driver2.Create(out_path, cols, rows, 1, gdal.GDT_Float32)

#Extent of template raster, which wilbe used for output raster
minx = gt[0]
maxx = gt[0] + cols*gt[1]
miny = gt[3] + rows*gt[5]
maxy = gt[3]

#Rasterize Polygons
gdal.Rasterize(output, flood_def, format=driver2, outputType=gdal.gdalconst.GDT_Float32,
               creationOptions=['COMPRESS=LZW', 'BIGTIFF=YES', 'blockxsize=256', 'blockysize=256'],
               outputSRS='EPSG:29902', outputBounds=[minx, miny, maxx, maxy],
               width=cols, height=rows, xRes=gt[1], yRes=gt[5], noData=ndv,
               targetAlignedPixels=True, allTouched='YES', burnValues=[1])

The error I am getting is:
---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-17-c5a7fba33323> in <module>
      1 #Rasterize Polygons
----> 2 gdal.Rasterize(output, flood_def, format=driver2, outputType=gdal.GDT_Float32,
      3                creationOptions=['COMPRESS=LZW', 'BIGTIFF=YES', 'blockxsize=256', 'blockysize=256'],
      4                outputSRS='EPSG:29902', outputBounds=[minx, miny, maxx, maxy],
      5                width=cols, height=rows, xRes=gt[1], yRes=gt[5], noData=ndv,

C:\ProgramData\Anaconda3\envs\env1\lib\site-packages\osgeo\gdal.py in Rasterize(destNameOrDestDS, srcDS, **kwargs)
   1161
   1162     if 'options' not in kwargs or isinstance(kwargs['options'], list) or _is_str_or_unicode(kwargs['options']):
-> 1163         (opts, callback, callback_data) = RasterizeOptions(**kwargs)
   1164     else:
   1165         (opts, callback, callback_data) = kwargs['options']

C:\ProgramData\Anaconda3\envs\env1\lib\site-packages\osgeo\gdal.py in RasterizeOptions(options, format, outputType, creationOptions, noData, initValues, outputBounds, outputSRS, transformerOptions, width, height, xRes, yRes, targetAlignedPixels, bands, inverse, allTouched, burnValues, attribute, useZ, layers, SQLStatement, SQLDialect, where, optim, callback, callback_data)
   1148             new_options += ['-optim', str(optim)]
   1149
-> 1150     return (GDALRasterizeOptions(new_options), callback, callback_data)
   1151
   1152 def Rasterize(destNameOrDestDS, srcDS, **kwargs):

C:\ProgramData\Anaconda3\envs\env1\lib\site-packages\osgeo\gdal.py in __init__(self, *args)
   3528     def __init__(self, *args):
   3529         """__init__(GDALRasterizeOptions self, char ** options) -> GDALRasterizeOptions"""
-> 3530         this = _gdal.new_GDALRasterizeOptions(*args)
   3531         try:
   3532             self.this.append(this)

TypeError: sequence must contain strings

Any idea where the problem is?


texts: ArcGIS Raster to Float tool ends up with different coordinate system

texts:
I am exporting a raster within a geodatabase to a flt file using Raster to Float tool in ArcGIS 10.7
The raster had an XY coordinate system of GDA 1994 MGA Zone 56.
When I checked the coordinate system of the Float file, it's GDA 1994 UTM Zone 56S.
Can someone explain why this is happening? I checked under "environment" to set the coordinate system to be "Same as Input" but this still happened.


texts: Import using wildcard in ArcGIS ModelBuilder?

texts:
I wish to take a shapefile export, which is kept in a folder location, and update a sqlite database. My problem is that the shapefile changes name slightly with each and every export.

Is there a way, in the ArcGIS ModelBuilder to handle importing the data with a wildcard?

The export shapefiles are the only ones stored in this folder - would be great if the add data handled Import *.shp.Edit;

Import shapefile, rename shapefile, make feature layer, create feature class based on schema of feature layer, export feature class to SQL database.


texts: Calculate cumulative distance between points in a dataset in ArcGIS 9.3 ArcInfo

texts:
I am working with several point datasets that I have received, with the points representing xy locations with elevation data for several levees in south Louisiana. I would like to be able to calculate the cumulative distance of each point from a point of origin.  By example, if my point of origin would represent 0, the distance to point 1 is 85 feet, the distance from point 1 to point 2 is 90 feet, the distance to from point 2 to point 3 is 70 feet, then the cumulative distance for point 1 would be 85 feet, point 2 would be 175 feet, point 3 would be 245 feet, etc.  And I would have to add the point of origin to several of the datasets for the purpose of consistency.

Thanks in advance for any assistance.


texts: Adding layer name and field entry where clause

texts: I am very new to the ArcGIS javascript API. I am creaing a web map using only one query layer in arcmap. If field 1 is not null I want to display x and if field 1 is null I want to display y. I did this previously by creating more than one layer in arcmap.However, I would like to only use this one query layer.


texts: XY points not showing up on map after adding new rows to attribute table

texts: In QGIS, I added two new rows, but the points are not showing up on my map even though I inputted their LAT and LONG's. The rest of the rows LAT/LONGs show up on the map (they were part of the original table I added). I already went into field calculator and specified that LAT was $y and LONG was $x...but that didn't solve the problem. How do I get the points to show up?


texts:
I am trying to use the Overpass Api with a python script to count the number of schools (and later other points of interest) for some area given as bounding box.
An example-query looks like this:

(node
['amenity'='school']
(48.58130844439305, 7.7423333727212675, 48.58435427306268, 7.7487015006498146);
way
['amenity'='school']
(48.58130844439305, 7.7423333727212675, 48.58435427306268,  7.7487015006498146);
rel
['amenity'='school']
(48.58130844439305, 7.7423333727212675, 48.58435427306268, 7.7487015006498146);
);
(._;>;);
out;

Using overpass-turbo, I verified that there are two schools in the specified area. What I need now is an automated way to extract that number. I found an example to get counts as output (here), but this only counts nodes/way/rels, and not "entities" I am interested in. The XML I get for above query contains 31 nodes and 1 way.
I thought about getting the base XML and count the number of "name"-tags as this seem to result in the number I need.

And here are my questions:

Is counting "name"-tags a sure way to get what I want, or will there be some examples where it won't work?
Is there any built-in way to get what I want? I would honestly be surprised if there is no faster/easier method.


texts: Create private access to qgis2web maps on Wordpress site

texts:
I want to use Wordpress site and qgis2web plugin to show maps to clients. According to my idea, a user logins to his private page and gets links to his private maps. But when I upload qgis2web folder to my site, everybody with the link can access this folder without authentication. And I can't give access to this folder to specific user, because this folder is not a wordpress page.

I can imagine two approaches of making qgis2web maps private:
1 - making qgis2web folder inaccessible by direct link and accessible from specific page. Does anybody know if it is possible?

2 - converting qgis2web folder (or at least index.html) to wordpress page. It is more than just copying content of index.html to blank page because you must replace all relative links with absolute links or copy content of all js and css. So does anybody know plugin (or third software) that can do it automatically?


texts: Shrinking coordinates of LAS file to fit in [0,1] using lidR

texts:
I normalised a LAS file's coordinates so that XYZ all range between 0 and 1.
Shift_Vox_X_ALS <- min(LAS_Vox_Test@data$X)
Shift_Vox_Y_ALS <- min(LAS_Vox_Test@data$Y)
Shift_Vox_Z_ALS <- min(LAS_Vox_Test@data$Z)

LAS_Vox_Test@data$X <- LAS_Vox_Test@data$X - Shift_Vox_X_ALS
LAS_Vox_Test@data$Y <- LAS_Vox_Test@data$Y - Shift_Vox_Y_ALS
LAS_Vox_Test@data$Z <- LAS_Vox_Test@data$Z - Shift_Vox_Z_ALS

Shift_Vox_N_X_ALS <- max(LAS_Vox_Test@data$X)
Shift_Vox_N_Y_ALS <- max(LAS_Vox_Test@data$Y)
Shift_Vox_N_Z_ALS <- max(LAS_Vox_Test@data$Z)

LAS_Vox_Test@data$X <- LAS_Vox_Test@data$X/Shift_Vox_N_X_ALS
LAS_Vox_Test@data$Y <- LAS_Vox_Test@data$Y/Shift_Vox_N_Y_ALS
LAS_Vox_Test@data$Z <- LAS_Vox_Test@data$Z/Shift_Vox_N_Z_ALS

After normalising the coordinates, although range(LAS_Vox_N@data$X) is:
[1] 0 1    # same for Y and Z

the LAS extent is still:
extent : 372179.5, 372194.5, 5828865.5, 5828880.5 (xmin, xmax, ymin, ymax)

But typing gives:
extent(LAS_Vox_N)
#> class      : Extent
#> xmin       : 0
#> xmax       : 1
#> ymin       : 0
#> ymax       : 1

I am not sure whether above behaviour is the reason why a saved las file (using writeLAS()) that is later opened (using readLAS()) produces erroneous y coordinates:
extent : 0, 1, 3681381.352, 3681381.352 (xmin, xmax, ymin, ymax)

writeLAS(LAS_Vox_Test, paste(Folder_AU_VOX_LAS_TYPE[AU],"/LAS_Vox_Test2_Norm.laz", sep=""))
LAS_OpenTest <- readLAS(paste(Folder_AU_VOX_LAS_TYPE[AU],"/LAS_Vox_Test2_Norm.laz", sep=""), select = "xyzp0")

Am I normalizing the coordinates  correctly so that all XYZ values are between 0 and 1? Obviously this removes the geographic position of the point cloud, which is what I need to do for this exercise. Is this allowed using LAS files in lidR?


texts: Leaflet L.geoJSON.getLatLng() is not a function

texts:
I am trying to get latlang of a marker so I can use it and fly to later.

map.flyTo(myLayer.getLatLng() ,16);

But why do I get getLatLng() not a function?

codepen example


texts: ISIS to GDAL to QGIS: Georeferencing Mars Images

texts:
I'm trying to use QGIS to view some images from the Mars Reconnaissance Orbiter Context Camera. I use ISIS3 to download the images and convert them to mapped .cub files, following the instructions here: http://isis.astrogeology.usgs.gov/IsisWorkshop/index.php/Working_with_Mars_Reconnaissance_Orbiter_CTX_Data

Now I'm trying to export these ISIS cube files to GTiff files, to view in QGIS, I can do this but the georeferencing doesn't seem to work. I've tried this in 3 different ways:

Use the ISIS command isis2std to generate tif (and tfw) files.

When I try to open the tif file in QGIS it asks me to choose a projection type (side question, is there a suitable projection for Mars? I always choose WGS 84, stupid I know but I don't know a more sensible option). Then once the raster is loaded it automatically displays as a grayscale image but my coordinates are not in degrees, they appear to be in pixels and have a range of 100s of thousands.

This method produces a tif image of ~100 Mb from an original .cub of close to 1 Gb.
Use gdal: gdal_translate -of GTiff input_mapped.cub output.tif.

This time I don't get asked to choose a projection type. The projection takes on 'USER:100002'. Hopefully this means that it is managing to read the appropriate Mars projection from the ISIS cub? The image displays as a solid grey box in grey scale this time, but I can see details by switching to a different colour mapping scheme. Again though my coordinates are not in degrees, they range in the hundreds of thousands again.

This method produces a GTiff with the same ~1 Gb size as the original .cub image.
Open the ISIS .cub file directly in QGIS.

Has basically the same effect as for method 2.

Because of the file sizes I'd prefer a variation of method 1, but with a suitable projection and georeferencing scheme for Mars.

Edit: Here is the georeference, and the output of gdalinfo on the original ISIS .cub file:

georeference (in python):

>>>import gdal
>>>from gdalconst import *
>>>fn = 'P01_001356_1747_XN_05S221W.map.cub'
>>>ds = gdal.Open(fn, GA_ReadOnly)
>>>gt = ds.GetGeoTransform()
>>>gt
(-22085.510544416, 5.1698292472885, 0.0, -234679.22885141, 0.0, -5.1698292472885)

gdalinfo (from terminal):

$gdalinfo file.cub
Driver: ISIS3/USGS Astrogeology ISIS cube (Version 3)
Files: P01_001356_1747_XN_05S221W.map.cub
       P01_001356_1747_XN_05S221W.map.cub.aux.xml
Size is 8526, 30111
Coordinate System is:
PROJCS["Sinusoidal Mars",
    GEOGCS["GCS_Mars",
        DATUM["D_Mars",
            SPHEROID["Mars",3396190,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Sinusoidal"],
    PARAMETER["longitude_of_center",138.45730511261],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-22085.510544415999902,-234679.228851410007337)
Pixel Size = (5.169829247288500,-5.169829247288500)
Corner Coordinates:
Upper Left  (  -22085.511, -234679.229) (138d 5' 1.74"E,  3d57'33.05"S)
Lower Left  (  -22085.511, -390347.957) (138d 4'56.04"E,  6d35' 7.46"S)
Upper Right (   21992.454, -234679.229) (138d49'45.19"E,  3d57'33.05"S)
Lower Right (   21992.454, -390347.957) (138d49'50.86"E,  6d35' 7.46"S)
Center      (     -46.528, -312513.593) (138d27'23.46"E,  5d16'20.26"S)
Band 1 Block=406x128 Type=Float32, ColorInterp=Undefined
  Min=0.033 Max=0.150
  Minimum=0.033, Maximum=0.150, Mean=0.080, StdDev=0.020
  NoData Value=-3.40282265508890445e+38
  Metadata:
    STATISTICS_MAXIMUM=0.15037615597248
    STATISTICS_MEAN=0.080181184414784
    STATISTICS_MINIMUM=0.033307410776615
    STATISTICS_STDDEV=0.020299974127104


texts: I have clipped a raster to exact polygon outline using image analysis (ArcMap 10). The resulting temporary clip layer is exactly how I want it to be. However, when I export/save the temporary clip layer it is no longer clipped exactly to the clipping polygon. Can anyone suggest why this might be and how to work around?


texts:
There were two posts from back in 2011 about the advantages of gvSIG and comparing QGIS and gvSIG. It's been 8 years since and I'm wondering if there might be an update since gvSIG and QGIS (at time of writing) are at v2.4 and v3.5, respectively.

Particularly, I'm interested to know what processes one software does better in than the other and how/why so (e.g. for Plugins, QGIS wins hands-down as the community is larger which has led to more plugins being developed).

Outside of the actual technology, I'm also curious as to how the community help has evolved. QGIS has a vast amount of help forums, websites, and tutorials, and I'm wondering if the response of the gvSIG community or if the resources available are equivalent or near equivalent to QGIS.


texts: Where is the GeoServer heatmap rendering transformation source code?

texts:
I'm trying to find the source code for the GeoServer Heatmap rendering transformation. I can't find anything relevant in the GeoServer GitHub repo.

This page in the developer guide includes a code snippet from the transformation but when I search for "public class HeatmapProcess implements GeoServerProcess" that page in the developer guide is the only result.

I'd like to learn how the transformation produces grid cell values as there's no documentation on the internal logic.


texts: Field calculator if and elseif statement with  values in ArcGis 10.1

texts:
Trying to figure out how to get this "Headwater" column to indicate "Yes" (or 1) if the value of "NextDown_Copy" is "" and "No" (or 0) if it has any other value (eg. 7100381010).

I attempted to try figuring out the Python code to do this using 1 and 0 coding (although I'd rather use Yes and No) (as illustrated below) but keep getting errors.
I consulted a few other posts (another one), but still wasn't able to find the problem.


texts: csv to feature classes not yet geospaced

texts:
Use case: preparing features for collector/fieldmaps for recording existing underground irrigation. Optimally my field tech would be able to go to a customer's house, see the customer is has[item named in table] and select it to add.
New to ArcGIS Pro. I've figured out how to add individual feature classes in a geodatabase and now am looking to import a CSV from our vendor in the format:

Item Name
Item Type
item cost

rotor head
point
$30

pipe
line
$.5 / ft

spray head
point
$60

etc
et cetera
maximae pecunia

The list is 1003 rows long so I need the import process rather than manual feature class addition, though if pressed I could write up an autohotkey to do that.
I've tried importing as a table, exporting the table as a dBase, and gdb, then importing them but they didn't show up when importing databases.
I've tried separating the points and lines into 2 CSVs, importing the point table, then export the point table as feature classes but I don't have coordinates making tool "XY Table to Point" not valid.


texts: Calculate Convex Hull based on attributes

texts:
I have a layer with lots of (spatially grouped) lines - lines belonging to a group share the same value in a field.
I would like create convex hulls around each of these groups using QGIS 3.8.
I learned that Vector Geometry->Convex hull is not attribute specific and the batch functionality works with a set of layers, not specific on attributes.

I tried the model builder, but I'm struggling - I would need to first create a set of objects grouped by that attribute (in SQL it would be something like "GROUP BY"), and then something like a "do convex hull for each in set".

Or create a number of new layers based on the attribute and use convex hull in batch mode? This would be ~100 new layers - could this be handled with QGIS?


texts:
I had to manually build M values for hundreds of levees (in the downstream direction) and wondering if there is a way to check that I didn't miss any levees. Is there a way I can look-up a table with all the "M" values other than clicking on on each individual levee?

Ideally I would like to be able to use Field Calculator to populate a new column with every segment's M value so I could see if there are any 0 values.


texts:
I am a GIS Intern for a city government, we are preparing to purchase ArcSDE and implement it to our current GIS setup. I'm curious as to what should be done before it is purchased. Our GIS consists of several file geodatabases for parcels, road centerlines, hydrology, etc. as well as several other purposes. As well as shapefiles for the previously listed.

Is there any sort of preprocessing/preparation work I should do before trying to implement ArcSDE?

I have tried looking around for tutorials and things and can't find any useful sources. I am asking about making changes to feature classes or to the file geodatabases themselves. We are getting Enterprise and plan on doing Check-In/Check-Out replication

We are currently using ArcGIS Desktop 10.2.1.


texts: How to properly assign spatial reference to converted AutoCAD data?

texts:
I know this is a frequent problem, I myself  confronted with this problem many times, yet somehow each Autocad dataset I received is different (although they're from the same area). To some data, the projection I gave them worked perfectly, and in this case it didn't.

So, I've recently received some AutoCad data. Converted them in ArcMap. They are missing spatial reference. But I know their location. They're suppose to be in a specific city in Romania.

So first thing I did, I gave them the Stereo 70 projection (Romania's national projection), I assumed this is the projection they were created in, the data gets somewhere in the middle of the Mediterranean Sea. Second thing I tried is is giving them WGS84, and they got somewhere in the Atlantic Ocean ;)

PS: I use ArcGIS 10.1

PS2: I've tried all 3 conversion methods (I know) in ArcMap.

PS3: The autocad data came with no spatial reference. I only know they should be in a specific city.


texts:
I have written a script to add a list of the aerial photo in sid format to a mxd. It works although the layers are displayed blue.

Is there a way to change the display to black and white with ArcPy?

[

import arcpy,os
    def addLayer(current_mxd, filePath):
        df = arcpy.mapping.ListDataFrames(current_mxd)[0]
        newlayer = arcpy.mapping.Layer(filePath)
        arcpy.mapping.AddLayer(df, newlayer,"BOTTOM")
    projectlocation = ""

    projectID = ""

    prep = "x" + "\n" + "x" + "\n" + "x"

    projectClientName = "HV" + "\n" + "Environmental Site Assessment"  + "\n" + ""

    #### add edmonton imagery to app folder #####
    filePathList = [r'V:\1102\active\110220065\gis\imagery\aep\1950_0135_041-ps.sid',r'V:\1102\active\110220065\gis\imagery\aep\1962_0818_077-ep.sid']
    coe = [r'r'Q:\Images\EDMONTON_IMAGERY\2012\EDMONTON_2012_Colour.sid']
    coeLayer = [r'Q:\Images\EDMONTON_IMAGERY\2014\2014.lyr',r'Q:\Images\EDMONTON_IMAGERY\2015\2015.lyr']

    current_mxd = arcpy.mapping.MapDocument(r'C:\Users\vlaw\Desktop\Application\template\mxd\historical_template.mxd')
    cwd =  os.getcwd()
    for frame in arcpy.mapping.ListDataFrames(current_mxd):
        frame.spatialReference = arcpy.SpatialReference(os.path.join(cwd, 'projection','NAD 1983 3TM 114.prj'))

        for lyr in arcpy.mapping.ListLayers(current_mxd,"",frame):
            if lyr.name == "Site_Boundary":
                lyr.replaceDataSource(r'V:\1102\active\110220065\gis\project_data',"SHAPEFILE_WORKSPACE","study_area")
                frame.panToExtent(lyr.getExtent())
            for a in filePathList:
                b = os.path.basename(a)
                d = b.replace(".sid","")
                year = b[0:4]
                photoID = b[5:13].replace('_',' ')
                addLayer(current_mxd,a)

                projectNumber = arcpy.mapping.ListLayoutElements(current_mxd, "TEXT_ELEMENT", "ProjectNumber or FigureID")[0]
                projectNumber.text = projectID
                projectLocation = arcpy.mapping.ListLayoutElements(current_mxd, "TEXT_ELEMENT", "Project Location Line")[0]
                projectLocation.text = projectlocation
                projectLocation.elementPositionY = 1.7813
                QAQC = arcpy.mapping.ListLayoutElements(current_mxd, "TEXT_ELEMENT", "Prepared By text block")[0]
                QAQC.text = prep
                clientProjectName = arcpy.mapping.ListLayoutElements(current_mxd, "TEXT_ELEMENT", "Client/Project Name")[0]
                clientProjectName.text = projectClientName
                clientProjectName.elementPositionY = 1.6322
                figureNumber = arcpy.mapping.ListLayoutElements(current_mxd, "TEXT_ELEMENT", "Figure Number1")[0]
                figureNumber.text = year
                reference = arcpy.mapping.ListLayoutElements(current_mxd, "TEXT_ELEMENT", "reference")[0]
                reference.text =reference.text.replace("REPLACE",photoID).replace("YEAR",year)
                mxdName = "fig_"+year+"_historical.mxd"
                ws = r'V:\1102\active\110220065\gis\figures\historical'
                pdf_folder= os.path.join(ws, "pdf")
                current_mxd.saveACopy(os.path.join(ws,mxdName))
                pdf_name = mxdName.replace(".mxd",".pdf")
                output = os.path.join(pdf_folder, pdf_name)
                arcpy.mapping.ExportToPDF(current_mxd, output, resolution =250)
                if lyr.name == d:
                    lyr.visible = False

                c = mxdName+ "\n" + "Exported map file: " + str(pdf_folder) + "\n"
                print c


texts: OS Translator II in QGIS failing to import - OGR version issues?

texts:
I'm trying to import mastermap into QGIS (on Mac OS) using the OS Translator II plugin. I'm a newbie!

I'm running QGIS 2.18.15, with GDAL/OGR 2.1.3, and the OS Translator II plugin was downloaded through the plugin manager.

QGIS was installed from http://www.kyngchaos.com/software/qgis, and I installed the four additional packages that were included in the bundle, including "1 GDAL Complete.pkg" - which appears to be /Library/Frameworks/GDAL.framework/Versions/1.11

Firstly when I run the OS Translator import it complained that ogr2ogr wasn't found. I'm not sure why it was not using the version installed with QGIS. I temporarily fixed it by modifying the import_task.py script to give it the explicit path:

joe:OSTranslatorII root# pwd
/Users/joe/.qgis2/python/plugins/OSTranslatorII
joe:OSTranslatorII root# diff -u import_task.py.orig import_task.py
--- import_task.py.orig 2018-02-11 11:39:06.000000000 +0000
+++ import_task.py  2018-02-11 08:23:56.000000000 +0000
@@ -34,7 +34,7 @@
                 cmd += arg + ' '
             print 'Starting %s' % cmd

-        self.process.start('ogr2ogr', self.args)
+        self.process.start('/Library/Frameworks/GDAL.framework//Versions/2.1/Programs/ogr2ogr', self.args)
     if not self.process.waitForStarted():
         raise Exception('Failed to start process. Please ensure you have gdal/ogr2ogr installed')

Then it could find the tool, but it complained that the ogr2ogr line options were wrong:

Args: ['-overwrite', '-lyr_transaction', '-f', 'PostgreSQL', '--config', 'PG_USE_COPY', 'YES', '--config', 'GML_GFS_TEMPLATE', '/var/folders/p4/rjqvjlbx2rz1cs9x74f0tb6m0000gp/T/1518284450.86.gfs', "PG:dbname='test' host='' port='5432' active_schema=os_mm_tmp user='postgres' password='XXXX'", '/vsigzip//Users/joe/.qgis2/python/plugins/OSTranslatorII/gfs/OS Mastermap Topography (v9) Pioneer.gz', '-lco', 'OVERWRITE=YES', '-lco', 'SPATIAL_INDEX=OFF', '-lco', 'PRECISION=NO', '--config', 'GML_EXPOSE_FID', 'NO']
Stdout: Usage: ogr2ogr [--help-general] [-skipfailures] [-append] [-update]
           [-select field_list] [-where restricted_where]
           [-progress] [-sql <sql statement>] [-dialect dialect]
           [-preserve_fid] [-fid FID]
           [-spat xmin ymin xmax ymax] [-geomfield field]
           [-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]
           [-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]
           dst_datasource_name src_datasource_name
           [-lco NAME=VALUE] [-nln name] [-nlt type] [-dim 2|3|layer_dim] [layer [layer ...]]

Note: ogr2ogr --long-usage for full help.

Stderr:
    FAILURE: Unknown option name '-lyr_transaction'

  Args: ['-append', '-lyr_transaction', '-f', 'PostgreSQL', '--config', 'PG_USE_COPY', 'YES', '--config', 'GML_GFS_TEMPLATE', '/var/folders/p4/rjqvjlbx2rz1cs9x74f0tb6m0000gp/T/1518284450.86.gfs', "PG:dbname='test' host='' port='5432' active_schema=os_mm_tmp user='postgres' password='XXXX'", '/tmp/2998191/2998191-summary.gml', '--config', 'GML_EXPOSE_FID', 'NO']
  Stdout: Usage: ogr2ogr [--help-general] [-skipfailures] [-append] [-update]
           [-select field_list] [-where restricted_where]
           [-progress] [-sql <sql statement>] [-dialect dialect]
           [-preserve_fid] [-fid FID]
           [-spat xmin ymin xmax ymax] [-geomfield field]
           [-a_srs srs_def] [-t_srs srs_def] [-s_srs srs_def]
           [-f format_name] [-overwrite] [[-dsco NAME=VALUE] ...]
           dst_datasource_name src_datasource_name
           [-lco NAME=VALUE] [-nln name] [-nlt type] [-dim 2|3|layer_dim] [layer [layer ...]]

I'm guessing that OS Translator II is implemented against a different version of ogr2ogr.

If anyone else has experienced this problem, is there a workaround?


texts: How to deal with missing values in Fishnet

texts:
I m working with ArcMap 10.2.2.

I have one dataset with red deer observation (Excel Sheet)
The data is joined to a polygon layer with administrative boudaries:
It looks like this

Blue are the red deer observations counted with zero (means data is available and zero animals are counted)
Yellow are areas where data didn't exist therefore I assigned the value -9999

I intersect my red deer data with a fishnet grid with cellsize 1000 -
No problem but when updating the fishnet - I summarize the red deer data by calculating the average.

I m getting wrong values in exactly this grid cells

Does anybody know how to deal with this problem.


texts: Matching two point clouds without loss (nearest neighbor, but more complicated)

texts:
I have two point clouds, one with demographic and one with electoral data. Points in both clouds represent the same empirical entity (say, a village, or an urban neighbourhood), but are not in the exact same lat/long (different datasources, GPS inaccuracies, etc). Moreover, one village can have several points in each cloud. I am looking for a cluster or nearest neighbor algorithm which allows me to merge both point clouds. I have thought about two strategies:

a) do a nearest neighbor matching, starting from the first point cloud and adding a variable to each point that directs me to the respective closest point in the second point cloud. The problem is that no method for nearest neighbor matching which I found (GRASS, SAGA, R) ensures that all points in the second cloud are picked at least once as this "respective closest point" for the first cloud. A neverending circle which could probably approximately be optimized - but I have no idea how. Any suggestions welcome.

b) the second strategy sees it as a clustering problem, in which I attempt to recreate the underlying empirical entity "village", which might contain multiple electoral and demographic points, but must contain at least one of each type. In other words, I would need a bottom-up hierarchical cluster algorithm which groups points together based on their spatial proximity, but simultaneously (!) ensures that each cluster contains at least one point from the electoral point cloud and one point from the demographic point cloud. Again: any suggestions welcome.

c) there are probably other solutions that I am unaware of - I am truly a beginner here...

It would be marvellous if somebody could point me in the right direction.


texts: How to export a shapefile through the Wikimapia API?

texts:
I am trying to download a kml file from Wikimapia containing information on quaters (districts), the very same information that is nicely available on Wikimapia. What I want to get at the end of the day is a shapefile containing as many polygons as there are districts in Milan to use in my GIS.

This is my URL for the query I pasted in my web browser.

http://api.wikimapia.org/?function=box&key=[mykey]&lon_min=8.936&lat_min=45.282&lon_max=9.421&lat_max=45.652&format=kml&category=quarters

From this I get a tiny kml.part file, that I then rename .kml. The file doesn't open in Google Earth or QGIS. If I open it with a text editor I can just few lines of xml, but no information concerning my geometries...

Am I doing something logically wrong?


texts: Leaflet Routing Machine Dynamically Generate Waypoints

texts:
I have a list of address that get geocoded. I want to then take those new lat/lon values and set them as waypoints in Leaflet routing machine so that an long route is generated. I have tried a loop which works, but it only does 2 points and then adds them to the map. The result is a lot of little routes on one map and not connected at all.  I cannot see a way to add dynamically add all waypoints in and then add the route to the map.  Is this even possible with Routing Machine??  Here is what I have tried:

HazardPoints is an array of lat/long values that I have cleaned up. I increment the loop by 6 each iteration (index 0 = first latitude value, index 1 = first longitude value, index 2 = hazard present[not used for this]) and then (index 3 = second latitude value, index 4 = second longitude value).. then loop jumps to the next set.

            function buildRoute(){
                var hazardPoints = finalPoints.toString().replace(/\n/g, ",").split(",");
                for(var i = 0; i < hazardPoints.length; i+=6){
                        L.Routing.control({waypoints: [L.latLng(hazardPoints[i], hazardPoints[i + 1]),L.latLng(hazardPoints[i+3], hazardPoints[i+4])],
                        routeWhileDragging: true}).addTo(map);
                }
                closeSideBar();
                    }


texts:
I have downloaded the Harmonized World Soil Data to get more information about global soil salinization and I want to process it in QGIS.

I have managed to open the raster layer (bil file) in QGIS and there is info in MDB files from the same database.

How can I connect the raster layer to the MDB files?


texts: Return empty dataframe when reading KML file of 'Sentinel 1 Acquisition Segments' using GeoPandas

texts:
I have downloaded the KML file of 'Sentinel 1 Acquisition Segments' from this url. I have written this code to read KML file into GeoDataFrame, but the code returns empty GeoDataFrame.
>>>import fiona
>>>import geopandas as gpd
>>> gpd.__version__
'0.8.1'
>>> fiona.__version__
'1.8.17'
>>>gpd.io.file.fiona.drvsupport.supported_drivers['KML'] = 'rw'
>>>fp = 'Sentinel-1A_MP_20201006T160000_20201026T180000.kml'
>>>gdf = gpd.read_file(fp, driver='KML')
>>>gdf
Empty GeoDataFrame
Columns: [Name, Description, geometry]
Index: []

Is there any code that I missed, or should I have to manipulate the KML file before?


texts: Looking for a fast, open-source raster cost-distance function to use in code

texts:
I've been scouring the internet for a fast, open-source cost-distance function that I can embed inside my code.  The functionality should basically work like the ArcGIS implementation where I can specify a source raster and cost raster and output a cost distance raster. That said, the implementation does not need to be as sophisticated as ArcGIS in that I don't need to necessarily input "rasters" with embedded geo metadata - a signature that takes simple numeric arrays will suffice since I can verify that the data overlaps correctly at the time I call the function.

My specific intent is to calculate the cost-distance to a single point in the center of a 1000 x 1000 cost raster based on a user click, so fast execution is very important.

A C# implementation would be ideal but I'll look at anything that's out there.

Does anyone know of an open-source library that supports this? Thanks for your help!


texts: Creating center-lines from amorphous-shaped polygons?

texts:
I would like to create road-center lines for amorphous-shaped images as in the example. void area is shown with white and solid area is shown with black. So I would like to create a center-line for the white area.

I have polygon version of the image also, but how can we produce center-lines for this image?

I've tried collapse dual lines to center-line command but it did not work.

When I tried to use dual lines, it does not create anything new or it just draws a very short line at the edge of the void area actually. What I want to do to create something similar to this.

I do not have ArcScan, I used ArcMap 10.3.


texts: Is it possible to delete only the graphical component

texts:
I use QGIS 2.10 and a Shapefile. This File contain a large attribute table and shows the data on the map like points. in short words,every thing is fine.

My question is, is it possible to delete only the graphical part not the whole database on qgis.
i would like to do that, maybe in future time i just want create new points based on new coordintes.

for example in MapInfo is it possible and in QGIS, too?


texts:
I'm using Leaflet to show some info for points:

var m1 = L.marker([20, 90], {icon: picture}).addTo(mymap);
m1.bindPopup('
Caption

Info: bla').openPopup();

How can you avoid that one popup is always opened when opening the HTML document?


texts: GIS newbie, How to georeference a jpeg metro map using OSS?

texts:
I have a metro map in jpeg format and I want to georeference it. My goal is to find GPS location of the metro stops and visualize them on google maps. Probably, the first step is to overlay the map over google maps area of interest to find the correct bounding box values. I was wondering, does anybody know a tool (web-based or standalone) for this.

Thanks,
-A


texts: Dead-ends are connected by osm2pgrouting

texts:
I've downloaded an extent of OSM data with QGIS's OpenStreetMap downloader, and want to use this data for with pgRouting. I can load it without problems as:

osm2pgrouting --f data.osm --conf mapconfig.xml --dbname db --username user --schema osm --clean

However, as you can see in the picture, the many of the dead-ends get connected to some seemingly random vertices. These dead-ends were created by cutting the OSM data with the extent polygon.

I'm using osm2pgrouting 2.3.2

Any tips on how to import OSM data cleanly?


texts: How to snap to intersection of lines with Classic snapping?

texts:
The question is as simple as the title.
When I am in Edit session with Classic Snapping enabled and want to snap to an intersection of lines, it does not work.
In order to achieve this, I have to disable the Classic Snapping, so ArcGIS can snap to the intersection.
Is there any way to snap to the intersection of lines with Classic Snapping enabled?


texts: Can Representations be used to Spatially Offset (Disperse) XY points for Web Mapping Application?

texts:
So I have a bit of an issue of points with the same x and y that I need to display separately from themselves. Here is some history of this issue:

I have a database of incidents that happen over a county. These incidents join to a parcel fabric with a one-to-many relationship. So several incidents occur on one parcel. I need the incidents to be shown as a point feature class for a web mapping application. So I added the x and y centroid of the associated parcel and created the point feature class. This resulted in several points being drawn on top of each other as shown in this picture:

I converted the symbology to a representation with a 10 meter minimum gap between points which resulted in exactly what I was looking for as show in this picture:

Now this works great in ArcMap, but when uploaded to my web mapping app the points can only be identified at there original x and y. So when you identify points that have been spaced you get no result, but when you identify the center of the parcel you get several results.

What I would love to happen is to somehow return the xy of these representations and change xy of each point to that. I've tried a few programming solutions, but I haven't been able to solve it yet.


texts: Working on databased-stored-geodatabase with a ArcGIS for Desktop standard license?

texts:
I'm looking to improve my GIS capabilities within ArcGIS for Desktop environment: I want to work with real geodatabase that is stored on a Microsoft SQL Server (SQLEXPRESS) and to have all abilities that are available when working on a File Geodatabase - that is: querying, editing, topology, relationship class and so on. Above that I want also to be able to use the Geodatabase Replication capabilities, i.e. working with a copy of the geodatabase and then syncing it in to the original geodatabase.

Is it possible to do so only having ArcGIS for Desktop Standard licence?

I've got a bit confused by ESRI documentation that mentions constantly that this ability is granted for this level of licencing but on the other hand requires a keycodes file that is granted when authorizing ArcGIS for Server Enterprise (!).

Standard licence example:

  Within ArcGIS for Desktop Standard and Advanced, the ArcCatalog application provides the ability for you to fully administer and manage enterprise geodatabases using SQL Server Express. This provides full enterprise geodatabase capabilities for up to a few users at a time. You set up and manage these enterprise geodatabases within ArcCatalog or the Catalog window. No extra software or database administration expertise is required.

Contradicting ArcGIS for Server example:

  To specify an Authorization File, browse to and choose the keycodes file that was created when you authorized ArcGIS for Server Enterprise.

In order to  achieve my goal I've did the following steps: (1) Installing Microsoft SQL Server 2012 SP1 (64-bit). (2) Installing the latest SQL Server Management Studio. (3) Installing the SQL Server Native Client downloaded from My Esri site.

I've succeed the following steps: (1) Connecting to the database server through ArcMAP's ArcCATLOG window. (2) Connecting to a pre-self-made database in the SQL Server. (3) Creating and importing feature classes. (4) Viewing and querying those feature classes.

I haven't succeed the following: (1) Creating new geodatabase within the database server nor using ArcMAPS's GUI nor using the GP tool. (2) Enabling geodatabase in this database. (3) Editing features within those feature classes.

It would be very helpful for me if you can resolve my main question and if the answer is true to point out what I've did wrong.


texts: Is there a way to create dot density maps with R package sf?

texts:
That is, is there a way to place dots at random positions inside a polygon with sf?

It's easy enough to do with maptools, using dotsinPolys, but since I like sf's "tidyness", I'd love to be able to make my dot density maps with it.


texts: Individual point attributes PostGIS Multipoint

texts:
I'm coming from a heavy Esri background and I'm trying to figure out how to create a multipoint feature in a PostGIS database that allows me to set a unique attribute for individual points in a multipoint feature. I've got data that is basically a huge number of XYZ points with an associated reading value (float) that are all grouped. I'd rather not create a database table with individual points as it will quickly baloon into the hundreds of millions of points. I've run databases with around 100 million rows before and it starts to get a little less performant past around 80 million rows in my experience.

I can always toss more hardware at it or implement table partitioning but I'd rather just use a multipoint and never run into issues to begin with.


texts: How to unlock imported GeoJSON vector layer?

texts:
I have a map on geojson.io, with some errors that I need to correct.

It has a feature to edit a polygon's node but it's too tedious and doesn't have "snapping" mode.  So I've imported my .geoJSON, .shp, and .topoJSON format files into QGIS.

And QGIS doesn't enable edit mode. Is it possible to do that or am I on the wrong path? If it's possible, how can I do it?


texts: Offsetting Advanced Digitizing lines from north in QGIS

texts:
In QGIS, I'm trying to draw some polygons (mostly rectangles) on a grid, with neat 90-degree right angles and such. All the polygons should be parallel or perpendicular to each other.
The problem is that the city I'm working on top of, Chicago, skews a little bit off a proper north-south line. For example, in the image below, Grant Park's field (yellow rectangle) is off just a bit from the red north-south line generated by Advanced Digitizing (and the red rectangle built from that line). Thus, anything I create off that red line won't match up with the park's grid.
Can I somehow "rotate" or "offset" the tool such that every line snaps to the yellow rectangle's lines, instead of compass north? If I add a bunch more red rectangles, their lines should all be parallel and perpendicular with the yellow one over Grant Park.
I tried rotating the entire QGIS canvas, but Advanced Digitizing still just points north-south. Or, if not through Advanced Digitizing, can I create my own guide/grid lines to snap other polygons to? Not really sure the proper way to do this.


texts: Tiled raster image to be placed in an OSM/ Leaflet.js map

texts:
I have a mapoverlay in a custom projection which I need to add as a layer to an Openstreet map within a website. Since I was not able to change the projection earlier, I decided to manipulate the raster image in Photoshop in order to make it "fit" onto the Openstreet base map.

Now I have a png file with the "correct" projection. My problem is that I now have to make tiles from this image which can be used as overlay and thus need to set the correct coordinates. I can find out the exact lat/lon values using OSM but how can I tell gdal2tiles where the top-left or bottom-right position of my image is?

My command line so far is:

$gdal2tiles.py -p 'raster' myimage.png

What I am missing is an idea on how to set coordinates. Do I have to use another tool prior to tile-making?


texts: How to change / switch mapsets in GRASS using cmd

texts: I am using GRASS 7.0.4 I have created a new mapset to the LOCATION using cmd in it. But unable to switch back to old mapset. can anyone guide me??


texts: How to change area inside mask to different RGB Value in ECW Format with FME

texts:
I have 3 Band RGB ECW Raster files and would like to change all the RGB values within the masked area to 0,0,0.

I am using the SpatialFilter to set the masking area and thought that I could use the RasterExpressionEvaluator(REV) to manipulate the rgb values. However I cannot find an expression which works.

I am using FME 2014 and setting the SpatialFilter to Mask = Filter, ECW File as Candidate with the Filter being "Intersects"and if does not pass, then the REV should work its magic. For each band in the REV (red8,green8,blue8) I set the the expression to

if ((A[0]+A[1]+A[2])>0, 0, A[0])

changing the last parameter for each band ...ie A1, A[2]...but this changes the whole image to 0,0,0. This doesn´t seem to be working.


texts: Shape looks like big Square in SQL Server 2008

texts:
I was provided with a polygon defined in WKT which has 4741 points and am able to insert it as a geography data type within a database in SQL Server 2008.  However, when I view it using the spatial viewer it takes up the entire viewer and looks like a large rectangle.  When I look at the same shape in Quantum GIS the shape looks like what I would expect (it is kind of a funky shape).  I executed the STNumGeometries() method on the geography shape in SQL Server and it returned 1 so SQL Server thinks it is valid, but does not appear to be able to correctly understand the shape.  Is there anything obvious that could be causing this issue?

The boundaries of the shape should be something like:

-89.25521,39.22638
-87.49851,38.84994
-87.97533,38.21419
-89.17156,38.11799
-89.25521,39.22638

What I'm seeing in the SQL Server Spatial viewer is:

-180, 81
180, 81
180, -90
-80, -90
-180, 81

Any assistance woudl be greatly appreciated.


texts:
I'm using gdaltransform to convert coordinates.

What's the best way to convert MGRS using gdaltransform given that it doesn't have an EPSG code?


texts: Scholarly articles that compare and contrast different thematic map types or discusses what thematic map types are better in certain situations?

texts: I'm writing a paper that involves some GIS work that I've done and as part of that I need to analyse an existing cluster map and explain why a heat map would be a better solution. I can provide plenty of reasons but they need to be backed up by some sort of literature review. Can anyone suggest a paper that would explain why certain map types would be better than others in different situations? I could use a web article in a pinch but it's better if it's some sort of research paper.


texts: Land cover change - need help with land cover classes

texts:
I'm very much new to using ArcGIS and to this forum. For my thesis, I have a raster map of Indonesia where I only look at three types of land cover retrieved from the ESA-CCI dataset. These involve: croplands, forest and grassland. I merged several land cover classes representing "croplands" as one value, the same I did for forest. In other words, I basically reclassified my land cover classes.

However, some land cover classes within "croplands" involve mosaic classes, where I had to assume my own cropland fraction (%). For example, for Mosaic cropland (>50%)/natural vegetation (tree,shrub,herbaceous cover)(<50%), I assume a cropland fraction of 60% and 40% respectively.

How do I set this fraction in my map?

Since the overall objective is to quantify how much forests has been converted (100%, one entire pixel basically) into croplands and grassland in a given year.


texts: Adding data to map and not saving it on local path using ArcGIS Pro?

texts:
Using ArcGIS Pro, I have loaded an APRX file that is located on a network share.

I added a layer to the map using Add Data.  The layer is named "USA_Detailed_Water_Bodies.lpk".

The added layer is stored locally in C:\Users\user\Documents\ArcGIS\Packages\USA_Detailed_Water_Bodies_41CAC1A1-63B9-4999-9149-7A6CB9942149\v104\dtl_wat.gdb

How can the layer be added to the map and stored with the map on the network share in order for other users to view the Water Bodies layer when the APRX file is loaded on another computer?


texts: Is it possible to run an arcpy execute SQL method without connecting to any database or server?

texts:
I want to simply process a layer by writing a tool script in a toolbox without having to connect any database, so I think I may have to set up a workspace environment in ArcMap, but I'm not sure if I need to do anything else.

Essentially I want to find a Python version of the ArcObjects IWorkspace.ExecuteSQL method.


texts: Unable to rename file geodatabase in ArcGIS running in a VM

texts:
I'm running ArcGIS 10.3 in a Windows 8.1 VM (using VMWare Player). I cannot rename file geodatabases in ArcCatalog. I make sure everything else is closed except for ArcCatalog, I create a New File Geodatabase, but when I try to name it anything else, it gives me the standard "General Function Failure".

I thought the problem was with the data being stored on a partition on my host which I access through a Windows mapped drive, but I've always been able to create/rename/delete folders, shapefiles, and everything else as normal, except for file geodatabases. I have full admin rights on my computer.

This led me to the conclusion that it has to do with the VM locking the gdb as soon as I create it, thereby preventing me from renaming it. However, even if I restart the VM, immediately open ArcCatalog and try to rename it, it still gives the same error.

I've tried renaming older file gdbs as well, with the same result. I've had this problem since running 10.2 on Windows 8. It's quite frustrating to have my neatly organised project folders filled with "New File Geodatabase.gdb", "New File Geodatabase_1.gdb" etc.

I can even rename feature classes and other objects within the file geodatabase, just not the actual file geodatabase. The question is: why cannot I rename file geodatabases in ArcGIS inside a Windows 8.1 VM using VMWare Player* on Windows 7?

*When I used this same VM for a time in VirtualBox, I did not have this issue.


texts: Getting subset with same cordinates for different rasters using gdal.translate?

texts:
I am a new into Python. I have processed 1000 NDVI rasters covering larger area. I need to burn those NDVI raster with Water mask. The water mask came for whole world I need to subset AOI and also want to subset the same AOI for the 1000 NDVI rasters. I want rows, column should be same so that I can burn with each other.

I have tried with https://automating-gis-processes.github.io/CSC18/lessons/L6/clipping-raster.html,
but it gives same rows column for all NDVI rasters but the no of columns and rows are changed when I subset the water mask with same cordinates.

I have tried with gdal.translate()

from osgeo import gdal

from shapely.geometry import box

minx, miny = 3.59, 41.14

maxx, maxy = 15.92, 49.27

bbox = box(minx, miny, maxx, maxy)

ds = gdal.Open('/scratch/test_reproject/20030416_1335_ndvi_reprojected_proj4.tif')
ds = gdal.Translate('/scratch/test_reproject/output/20030416_1335_ndvi_reprojected_proj4_crop_gdal.tif', ds, projWin = [3.59, 41.14,  15.92, 49.27])
ds = None

from osgeo import gdal

from shapely.geometry import box


texts: GeoEvent Processor for Server Service Stopping

texts:
I have recently began using the GeoEvent Processor for server, It will run for one minute and then stop. I can't really intepret the logs to see what the problem is, and I am hoping a better eye than mine can help me out. The log is from when I started the service to when it failed.

http://pastebin.com/M63dWcDW


texts: Read numerical value of variable in NetCDF4 using Python

texts:
I need to undestand how I can get the numerical value of one channel IR_108 ( variable) in NetCDF4 and some other data:

import numpy as np
from netCDF4 import Dataset

 path = '/home/data/Mmultic3kmNC4_msg04_201905080200.
    ...: nc'

In [78]: f=netCDF4.Dataset(path,'r')

 f.variables.keys()
Out[80]: odict_keys(['time', 'dtime', 'IR_016', 'commentaires', 'satellite', 'geos', 'ImageNavigation', 'GeosCoordinateSystem', 'Y', 'X', 'Albedo_to_Native_count_IR_016', 'IR_039', 'Temp_to_Native_count_IR_039', 'IR_087', 'Temp_to_Native_count_IR_087', 'IR_097', 'Temp_to_Native_count_IR_097', 'IR_108', 'Temp_to_Native_count_IR_108', 'IR_120', 'Temp_to_Native_count_IR_120', 'IR_134', 'Temp_to_Native_count_IR_134', 'VIS006', 'Albedo_to_Native_count_VIS006', 'VIS008', 'Albedo_to_Native_count_VIS008', 'WV_062', 'Temp_to_Native_count_WV_062', 'WV_073', 'Temp_to_Native_count_WV_073'])

In [81]: IR=f.variables['IR_108']

In [82]: IR
Out[82]:
<class 'netCDF4._netCDF4.Variable'>
int16 IR_108(ny, nx)
    standard_name: toa_brightness_temperature
    long_name: 10.8 microns infra-red channel MSG Imager
    units: K
    scale_factor: 0.01
    add_offset: 273.15
    valid_min: -27315
    valid_max: 32767
    _FillValue: -32768
    available: yes
    comment: channel 9
    nuc: 931.122
    alpha: 0.9983
    beta: 0.6256
    CalibrationParametersOrigin: rttov
    _CoordinateSystems: GeosCoordinateSystem
    grid_mapping: GeosCoordinateSystem
    resolution: 3.0
    Offset : radiance = offset + CN x slope: -10.456819486590666
    Slope : radiance = offset + CN x slope: 0.2050356762076601
unlimited dimensions:
current shape = (3712, 3712)
filling on

print(IR.dimensions)
('ny', 'nx')

fisrt question The dimensions (nx,ny) is the row and column? or the numerical value corresponds to coordinates x;y ?

What I would like to understand from the information on the variable IR_108 how I can have the numerical values in this channel and also the geographical coordinates X Y or lat and lon.

The idea that I have points x, y ( 2150,525) which corresponds to the line and column and I will have to extract from these x y the numerical value of the variable in this pixel. For example, the pixel 2150 row 525 corresponds to the value IR_108 = 4.7.


texts: How to express not within/inside predicate in spatial database with good performance

texts:
In our application we want to express some spatial predicates as objects not within some areas, but unfortunately it has very bad performance. The following SQL is generated by our application:

SELECT *
FROM geom_table
WHERE NOT (SDO_INSIDE(geom_column, SDO_GEOMETRY(...), srid) = 'True'
AND geom_column IS NOT NULL);

I've been searching on internet for a while but surprisingly fail to find any information talking about it. I think it's not a non-reality use case, or is it? I wonder if there is a better way to express the SQL statement?


texts: Resolve Error 000210 when using Scratch Workspace?

texts:
Now I'm trying to make Python script to resolve some problem.

I need too place all out features in Scratch Workspace (Later this script will be used in web GIS and intermediate features mustn't be saved).
But when I'm trying to do this i have and 000210 error.

How can i resolve this?


texts:
I just wrote a code for an online application through GEE. The application works when I click on the "run" button within the GEE interface.  However, there is no way to run the code in the app interface because there is no "run" button.
Is there any way to add a button to run the code in the app interface?


texts: GIS grid linear interpolation for a given point

texts:
Suppose that I have a GIS grid.

I need to do linear interpolation for a given point located inside a cell of the GIS grid.

The locations and values of the vertices of each cell in the grid are known.

I am also given the location and value of the centroid of each cell.

For a given point, I need to use linear interpolation to find the value.

If the given point is located on the centroid, I do not need to do any linear interpolation. I just pick up the value of the centroid.

But, if the given point is close to the centroid but not located on it, I want the value of the centroid to play an important role in the interpolation with the 4 vertices of the cell and the centroid.

How to implement this idea ?

I do not want to define a fixed value as a radius. If the given point distance from the centroid is less than the centroid, the centroid value will be used as the given point value. This may cause some "jumps" in the interpolation.


texts:
Following the OpenLayers tutorial I've got a working custom control with working tooltips for the default Controls as shown in this tutorial. Is it possible to add a tooltip to one or more custom controls in a simple way or to inherit that behaviour? Regarding the API, the custom control doesn't seem to inherit any tooltip behavior at all, yet the subclasses all have default Tooltip texts. If I, however, put in the same parameter for the tooltip as used in a default control, nothing happens.

// namespace
window.app = {};
var app = window.app;

app.customControl = function (opt_options) {
    var options = opt_options || {};

    var button = document.createElement("button");
    button.innerHTML = "";

    var custCtrlFct = function (e) {
        // Function logic
    };

    button.addEventListener('click', cstCtrlFct, false);

    var element = document.createElement('div');
    element.className = "custom-button ol-unselectable ol-control";
    element.appendChild(button);

    ol.control.Control.call(this, {
        element: element,
        target: options.target
    })
};
ol.inherits(app.customControl, ol.control.Control);

The custom control gets added to the map at creation.

var map = new ol.Map({
    controls: ol.control.defaults({
        attributionOptions: /** @type {olx.control.AttributionOptions} */ ({
            collapsible: false
        })
    }).extend([new app.customControl({tipLabel: "Tooltip"})]),
//[...]
});

Do I have to re-implement the whole Tooltip behavior by myself or is there a way to get OpenLayers to use the Tooltip behaviour of the built-in controls?


texts:
How to increase GeoWebCache speed, Am using geoserver 2.7.2 for tile cache and data store creating from a SQL Server database it will take 1 day and above to complete cache and it terminated several time,am tring this using postgreSQL DB it will comparative faster than SQL server and it never terminate, what the reasons for that?

How increase GeoWebCache speed using data store from a SQL Server database?


texts: Rasterize right method for true value of point grid?

texts: Is rasterize an interpolation or does it just burn in the "true" value of the point? I have a grided data that is evenly spaced and dense. I would like to just have a raster with the "real" value of the point grid and where there is no value it should be Nodata. Is this the method to go if I only whant the values of the point?


texts: Why are values coming in as  when importing table from Excel?

texts: I'm importing a CSV file into ArcGIS Pro, I've used geoprocessing tool XY Table to Point but values have come in null. Column type is double in ArcGIS Pro. Highlighted is the column which imports null from Excel.


texts:
I need to know the minimum hardware requirements to run GRASS 7, at the request of a Journal editor. There is plenty of information in the Wiki about software requirements, but none on hardware.

All I could find is this unsigned message from 20 years ago at GisLounge:

  8 Mbytes of memory (of course, more is better..)

  100 Mbytes of free disk space

  ~40 mb for executables,

  ~40 mb for source code (which you can ignore if you merely install the
  Linux binaries)

  ~? for data (the veritable bottomless pit can be filled with data, if
  you so choose)

Naturally, things might be a little different by now.


texts: Update QGIS 2.6 with a bug fix (revision 11c2bae7)

texts:
How do I use a cpp file to update my QGIS?

I'm talking in this case of the issue in this link.

Can anyone help?


texts:
I am always in doubt as to where to categorize GIS. For instance, when you use a search engine to look for a job, you are asked to choose the category your job falls in. Or when you search for a book on an online library:

The above list contains the common major fields you are asked to choose from.

What's the most relevant field to categorize GIS?


texts: Lock the nativity of the QGIS using PyQGIS

texts:
How to lock the nativity of the QGIS using PyQGIS

icon = QtGui.QIcon(folPath + '\Dsmsoft.png')

pixmap = icon.pixmap(302, 302, QtGui.QIcon.Normal, QtGui.QIcon.On)
icon = QtGui.QIcon()
icon.addPixmap(pixmap, QtGui.QIcon.Normal, QtGui.QIcon.Off)
self.dlg.setWindowIcon(icon)


texts: OpenLayers WFS implementation in java for editing and saving data

texts: I am  new user in openlayers and want to implement WFS layer to edit and save data to postgis tables.am using netbeans as development environment and m very new in that.Can anyone to help me regarding to this?It would be good if any running example found.


texts: Change point symbol to an existing single marker symbol in PyQGIS

texts:
I have a code that takes a polygon layers and returns a temp point layer, I want to change the symbology for the point layer to the "topo pop capital" point symbol through Python.

fn = "C:/Users/Public/Desktop/OSGeo4W/bigBoi.shp"

layer = QgsVectorLayer (fn)

feats = layer.getFeatures()

vl = QgsVectorLayer("Point", "temp", "memory")
pr = vl.dataProvider()

for feat in feats:
    geom =feat.geometry()
    geomSingleType = QgsWkbTypes.isSingleType(geom.wkbType())
    if geom.type() == QgsWkbTypes.PolygonGeometry:
        if geomSingleType :
            x = geom.asPolygon()
            for i in x:
                for y in i:
                    for q in y:
                        f = QgsFeature()
                        f.setGeometry(QgsGeometry.fromPointXY(q))
                        pr.addFeature(f)
                        vl.updateExtents()
                        QgsProject.instance().addMapLayer(vl)

        else:
            x = geom.asMultiPolygon()
            for i in x:
                for y in i:
                    for q in y:
                        print(q)
                        f = QgsFeature()
                        f.setGeometry(QgsGeometry.fromPointXY(q))
                        pr.addFeature(f)
                        vl.updateExtents()
                        QgsProject.instance().addMapLayer(vl)

I found the setSymbol fucntion that takes a symbol and renders it on a layer like so:

symbol = QgsMarkerSymbolV2.createSimple({'name': 'square', 'color': 'red'})
layer.rendererV2().setSymbol(symbol)

I also found that you can "copy symbol" by right clicking on an existing layer with the desired symbol, but that returns an HTML-like code which for some reason I can't paste here.

Is there a way to create the symbol with the copy symbol output as an object?
Is the topo pop capital a symbol layer? And if so how do I appendSymbolLayer() if I don't know how to reference it.
Can I call to the symbols from the single symbol favorites?
Or maybe export the symbol to get it in Python template?


texts: Reprojecting WGS 1984 Web Mercator (EPSG:3857) in Python with GDAL

texts:
I am reprojecting rasters in python using GDAL. I need to project several tiffs from geographic WGS 84 coordinates to WGS 1984 Web Mercator (Auxiliary Sphere), in order to use them later in Openlayers together with OpenStreetMap and maybe Google maps. I am using Python 2.7.5 and GDAL 1.10.1 from here, and transforming coordinates using advices from here (my code is below). In short, I imported osgeo.osr and used ImportFromEPSG(code) and CoordinateTransformation(from,to).

I tried firstly EPSG(32629) which is UTM zone 29 and got this projected raster (more or less fine), so the code seems to be correct:

Then I used EPSG(3857) because I've read this and this questions and found that it is the correct recent valid code. But the raster is created with no spatial reference at all. It is far away up in WGS 84 data frame (but will be ok if I switch the data frame to Web Mercator).

With EPSG(900913) the output is georeferenced, but shifted about 3 raster cells to the north:

When I reproject the raster using ArcGIS (export in WGS_1984_Web_Mercator_Auxiliary_Sphere) the result is nearly fine:

And when I use the old code 102113 (41001,54004) the result is perfect:

The summary of my tests using all codes:

3857: far away up (missing georeference)
3785: far away up (like 3857)
3587: far away right
900913: slightly jumped up
102100: python error
102113: perfect
41001: perfect
54004: perfect
ArcGIS (web merc. aux.): good

So my questions are:

Why does the correct EPSG code give me wrong results?
And why the old codes work fine, aren't they deprecated?
Maybe my GDAL version is not good or I have errors in my python code?

The code:

    yres = round(lons[1]-lons[0], 4)  # pixel size, degrees
    xres = round(lats[1]-lats[0], 4)
    ysize = len(lats)-1  # number of pixels
    xsize = len(lons)-1
    ulx = round(lons[0], 4)
    uly = round(lats[-1], 4)  # last
    driver = gdal.GetDriverByName(fileformat)
    ds = driver.Create(filename, xsize, ysize, 2, gdal.GDT_Float32)  # 2 bands
    #--- Geographic ---
    srs = osr.SpatialReference()
    srs.ImportFromEPSG(4326)  # Geographic lat/lon WGS 84
    ds.SetProjection(srs.ExportToWkt())
    gt = [ulx, xres, 0, uly, 0, -yres]  # the affine transformation coeffs (ulx, pixel, angle(skew), uly, angle, -pixel)
    ds.SetGeoTransform(gt)  # coords of top left corner of top left pixel (w-file - center of the pixel!)
    outband = ds.GetRasterBand(1)
    outband.WriteArray(data)
    outband2 = ds.GetRasterBand(2)
    outband2.WriteArray(data3)
    #--- REPROJECTION ---
    utm29 = osr.SpatialReference()
#    utm29.ImportFromEPSG(32629)  # utm 29
    utm29.ImportFromEPSG(900913)  # web mercator 3857
    wgs84 = osr.SpatialReference()
    wgs84.ImportFromEPSG(4326)
    tx = osr.CoordinateTransformation(wgs84,utm29)
    # Get the Geotransform vector
    # Work out the boundaries of the new dataset in the target projection
    (ulx29, uly29, ulz29) = tx.TransformPoint(ulx, uly)  # corner coords in utm meters
    (lrx29, lry29, lrz29) = tx.TransformPoint(ulx + xres*xsize, uly - yres*ysize )
    filenameutm = filename[0:-4] + '_web.tif'
    dest = driver.Create(filenameutm, xsize, ysize, 2, gdal.GDT_Float32)
    xres29 = round((lrx29 - ulx29)/xsize, 2) # pixel size, utm meters
    yres29 = abs(round((lry29 - uly29)/ysize, 2))
    new_gt = [ulx29, xres29, 0, uly29, 0, -yres29]
    dest.SetGeoTransform(new_gt)
    dest.SetProjection(utm29.ExportToWkt())
    gdal.ReprojectImage(ds, dest, wgs84.ExportToWkt(), utm29.ExportToWkt(), gdal.GRA_Bilinear)
    dest.GetRasterBand(1).SetNoDataValue(0.0)
    dest.GetRasterBand(2).SetNoDataValue(0.0)
    dest = None  # Flush the dataset to the disk
    ds = None  # only after the reprojected!
    print 'Image Created'


texts: Problem printing maps to PDF with MapInfo

texts:
I'm having trouble printing maps to PDF using MapInfo Professional 12.0.1 at work.

The error on the PDF is:

I am using File -> Print to PDF -> MapInfo PDF Printer Version 12.0. The PDF needs to be a layered PDF due to the amount of information, but also add the attribute data - I believe these two options are causing the problem? Maps have exported to PDF in the past when these options are not ticked.

Regarding the print settings, I have enabled the options I found online which are meant to fix the issue (enabled postscript simulation, multi-language support, font embedding, tried multiple compressions to reduce image size including 72DPI)

Is there something I am missing to get this to work? What is annoying is the exact same settings worked on a trial export at A3, but not now (only added a legend in).


texts: Merge/join only the lines which are snapped to vertices in QGIS

texts:
I have polygons drawn with AutoCAD as Mline (Multiline type) which are rendered in QGIS as 8 segments for each feature. An example of one selected feature is shown in the picture below:

I want to clean this layer thus to only have the connected segments as polylines. An illustration of the expected result is shown below:

It is not possible to convert them with AutoCAD.

I used the "Multipart to Singleparts" processing to obtain 8 features instead of 1 feature composed of 8 separated segments. Here is my problem, I do not find the process to keep only the connected segments and remove the isolated ones.

Then I tried several procedures:

play with the "v.clean" tool from GRASS, without success.
use the "Shortest path" processing, which is not a solution, as I have too many features.
"reconstruct" polylines by generating some points along lines and connect them automatically by creating new lines lines, with the help of some plugins.

Finally, to circumvent this issue, I also tried to apply the "Buffer" processing on the 8 segments to obtain overlapping polygons. Then, the "Dissolve" processing was applied on the overlapping polygons in order to obtain one polygon describing the former lines originally drawn as multilines. Finally, the use of the "Skeleton/Medial Axis" function of the HCMGIS plugin led to a result close to the one I expected. Unfortunately, three major problems are encountered in following this procedure: loss of accuracy, unexepected small polylines generated at each curvature and a very heavy procedure (crash) for my computer (more than 3000 segments).

Does someone have an idea?


texts: GIN/GiST Full Text searches through openlayers, geoserver and postgres

texts:
I have a a web app that automatically maps recent tweets and allows you to search for certain keywords. it uses Postgres for the database, geoserver as the server, and openlayers as the cartographic library. Right now tweet keywords are done using the ILIKE query.

This works pretty well if the keyword i'm looking for is a commonly used word, like 'love' , 'friend', 'OMG', etc. but less common words like "geospatial" need to search through a lot more data in order to find the last n instances of the word. This is dead slow.

To combat my slowness I'm building a GIN index on my tweet field in postgres. Two days later and i'm still waiting for the index to complete. Once it's built though, I'll try a few queries using SQL view parameters and hopefully this will dramatically speed things up.

Even if this approach works though, it's not going to be all that useful if the index takes so long to create -- the maintenance of it will not be able to keep up will the incoming flood of data. I think I will try out the GiST index next and see how long it takes as I understand it is much faster to this build index.

Beyond this, what can I do next? Do options like Solr work with Geoserver?? Are there any common approaches to this type of problem and data size/rate?


texts: Find/Replace Double Quotes using Python Parser of ArcGIS Field Calculator?

texts:
I'm trying to use the Field Calculator (ArcMap 10.0) to strip out double quotes in an address field (I didn't put them in there, believe me).  I've tried many incantations but am still unable to even get off the ground with a calculation. yes, both fields are string, and both are of acceptable length.

This calculation:

Results in this:

With this error message in the Results dialog:

And I repeat, I didn't put them in there.


texts: Dynamically change GeoServer SLD style on WMS call using ENV

texts:
I need to chage layer style depending on attribute value on wms call. Like this

<Fill>
  <CssParameter name="fill">
    <ogc:Function name="Recode">
     <!-- Value to transform -->
     <ogc:Function name="strTrim">
       <ogc:PropertyName>ward_no</ogc:PropertyName>
     </ogc:Function>

     <!-- Map of input to output values -->
     <ogc:Literal>28</ogc:Literal>
     <ogc:Literal>#6495ED</ogc:Literal>

     <ogc:Literal>10</ogc:Literal>
     <ogc:Literal>#B0C4DE</ogc:Literal>

     <ogc:Literal>6</ogc:Literal>
     <ogc:Literal>#F5DEB3</ogc:Literal>

     <ogc:Literal>11</ogc:Literal>
     <ogc:Literal>#F4A460</ogc:Literal>

     <ogc:Literal>14</ogc:Literal>
     <ogc:Literal>#87CEEB</ogc:Literal>

   </ogc:Function>
 </CssParameter>
</Fill>

But the attribute value [28,10,6,11,14] cal also be change dynamically. By using ENV value I can change color of each and every attribute, like this

<ogc:Function name="property">
    <ogc:Function name="env">
         <ogc:Literal>28</ogc:Literal>
         <ogc:Literal>#6495ED</ogc:Literal>
    </ogc:Function>
 </ogc:Function>

But how to change all input attribute?
Or is their any way to create a attribute based style on wms call??

The attribute value list can also increase or decrease also depending on request


texts:
I have built a tool that allows my users to perform geocode requests against the World Geocoding Service provided from Esri.

I would like to extend my tool to have an autocomplete feature, similar to how this can be achieved with the Autocomplete class from the Google Geocode API.

Although there is not a similar class available at the Esri service end, I was thinking that you could still hit the geocode service with a GET request, after listening for a slight delay after each keystroke from the end-user in something like a combo-box, and just return the top 5 results based on score.

I am looking for an answer that provides some code (can be pseudocode, or any language) on how to accomplish this, or if this workflow is even achievable.  Something similar to this, but against the Esri Geocode service as opposed to map services.


texts: Adding ECW support in GeoServer OpenGeo Suite?

texts:
I have installed opengeo suite 4.8 with GeoServer:

" 2.8-SNAPSHOT

Git Revision5ba955074e0c5b78a5c9576dba955bec93008889

Build Date14-Jan-2016 20:47

GeoTools Version14-SNAPSHOT (rev 109456cbb2cf31ce318330bad5d0176e3d04e4a9)

GeoWebCache Version1.8-SNAPSHOT (rev 1.8.x/854b89c13d79df5a232c32245f956cfb70de56f5) "

I am trying to add ECW support to this version of geoserver but no luck; don't know where can i download the right plugin, and don't know if that is possible on this version of geoserver. I would like to add other plugins if its possible in this version.


texts:
I would like to make improvements in the resolution of satellite images using pansharpening. After found this website, which explains how Planet Labs do it.

For a University project I would like to try different algorithms for pansharpening, so I could follow that guide but not using the script gdal_pansharpen.py. But here I have a problem: I couldn't find the script gdal_pansharpen.py so I don't know what is written in there.

I would choose Python for programming the different algorithms (IHS, BT, ...) and see which one is better.


texts:
I'm trying to write a script that looks into a dataset and list all the elements in the dataset. The dataset comprises of Feature classes and FGDB Toppology. When I uses ListFeatureClasses, it outputs all the feature classes in the dataset,so I know the code works, but I want to list also the topology file. ListFiles does not work. Which of the List functions should I use to list topology?
This is the code I have so far:

import arcpy
import os
PATH1 = r"C:\Users\Daimon Nurse\Desktop\DFMPROJECT\DFM.gdb"
PATH2 = r"C:\Users\Daimon Nurse\Desktop\DFMPROJECT\Complete Final Capstone Map\Z2 Map File Geodatabase.gdb"
arcpy.env.workspace = PATH2
arcpy.env.overwriteOutput=True
for data in arcpy.ListDatasets():
    if data == "AllZones":
        print data + "1"
        path = os.path.join(PATH2,data)
        print path
        arcpy.env.workspace = path
        for file in arcpy.ListFiles():
            print file

If I change arcpy.ListFiles() to arcpy.ListFeatureClasses(), the code works as I stated above.


texts: Layering WMS or ArcGIS Rest on MapBox IOS SDK

texts:
We have some customer's who have facility data they are publishing through their own WMS or ArcGIS/REST services.

We're exploring the idea of building an IOS app that can display their facilities using MapBox's IOS SDK.

We've had no trouble configuring a set of landbase layers in MapBox, and displaying them in our demo app. We've not yet tried uploading sample facility data onto the MapBox site, but we don't expect to have any difficulty with that.

We expect some of our customers will want to upload their facility data into MapBox, so that they can take advantage of MapBox's offline mapping features. But a good many of them won't.

For those that don't, we're trying to figure out how to access our customers' facility services directly.

MapBox IOS SDK supports custom raster styles: Custom raster style

The idea is you create a custom style in json, and in it you specify a template for the tile URLs.

{
  "version": 8,
  "name": "Raster Tiles",
  "sources": {
    "yourTileLayer": {
      "type": "raster",
      "tiles": [
        "https://1.tile.server/{z}/{x}/{y}.png",
        "https://2.tile.server/{z}/{x}/{y}.png"
      ],
      "tileSize": 256
    }
  },
  "layers": [{
    "id": "yourTileLayer",
    "type": "raster",
    "source": "yourTileLayer",
    "paint": {
      "raster-fade-duration": 100
    }
  }]
}

Those templates look suspiciously like the templates use by Leaflet's TileLayer: TileLayer

L.tileLayer('http://{s}.tile.osm.org/{z}/{x}/{y}.png?{foo}', {foo: 'bar'}).addTo(map);

Leaflet has a TileLayer.WMS class, that derives from TileLayer, that provides an interface that's usually used for displaying WMS on a Leaflet map.

Is there a template I could pass to Leaflet's TileLayer class, that would retrieve tiles from a WMS server, without using Leaflet's TileLayer.WMS class?

Can anyone tell me how to build it?

Similarly, ESRI has a Leaflet plugin to display ArcGIS/REST data on a Leaflet map. Is there a template I could pass to Leaflet's TileLayer class, that would retrieve tiles from an ArcGIS/REST server, without using ESRI's plugin?

Can anyone tell me how to build that?

My hope is whatever URL I might build that works with Leaflet's TileLayer class will also work with MapBox's IOS SDK.


texts: Customize Toolbar - customize location of savings

texts:
In Where does QGIS Plugin "Customize Toolbar" save its changes? it was pointed out, that "Customize Toolbars" saves the changes in a file called ".CustomToolBars". This file is stored in the home directory of the user, so the path is something like "c:\users\xyz"

I am running a mobile version of QGIS and I would like to change the location  of this file to the current directory where the plugin is installed (or any other directory  on the USB drive). I use this mobile version for teaching and need to have the path set dynamically so the file will be found no matter  of the path of the mobile QGIS.

I have only very little knowledge of programming (and especially Python), but I had a look at the python code and I think that the code lines for defining the location are to be found in line 50 and 51 of "CustomToolbar.py":

[50] self.userhome = os.path.expanduser('~')
[51] self.filepath = self.userhome + '//.CustomToolBars'

In order to set the location dynamically to the current directory of the plugin I changed it as follows:

[50] self.userhome = os.path.dirname(__file__)

But when starting QGIS it fails to load the plugin and gives the following errors:

Traceback (most recent call last):
  File "H:/GIS/usbgis/apps/qgis2/apps/qgis-ltr/./python\qgis\utils.py", line 219, in startPlugin
    plugins[packageName] = package.classFactory(iface)

File "H:\GIS\usbgis\config_lite/python/plugins\CustomToolBar\__init__.py", line 10, in classFactory
    from CustomToolbar import CustomToolbar

  File "H:/GIS/usbgis/apps/qgis2/apps/qgis-ltr/./python\qgis\utils.py", line 478, in _import
    mod = _builtin_import(name, globals, locals, fromlist, level)

  File "H:\GIS\usbgis\config_lite/python/plugins\CustomToolBar\CustomToolbar.py", line 52
    self.userhome = os.path.dirname(__file__)
   ^
IndentationError: unexpected indent

Actually, I am a bit lost and don't quite know how proceed further.


texts:
We are now using QGIS 3.12.
Our basic python script from earlier QGIS version  to export to Mapinfo and mid/mif seems to no longer work.
I understand QGIS 3+ has upgraded this functionality. and i have no idea where to start to fix this.

Can someone please advise how to run the below script, or advise how to upgrade it also?
Can someone please review?

from qgis.gui import*
from qgis.core import *
import os

def run_script(iface):
    PATH = QgsProject.instance().readPath("./")
    MifPath = PATH + '/Submit/MID_MIF/'
    ShapePath = PATH + '/Submit/ESRI_SHAPE/'

    if not os.path.exists(MifPath):
        os.mkdir(MifPath)
    if not os.path.exists(ShapePath):
        os.mkdir(ShapePath)

    for vLayer in iface.mapCanvas().layers():
        QgsVectorFileWriter.writeAsVectorFormat( vLayer,
        MifPath + vLayer.name() + ".mif", "utf-8",
        vLayer.crs(), 'Mapinfo File', datasourceOptions='FORMAT=MIF')

    for vLayer in iface.mapCanvas().layers():
        QgsVectorFileWriter.writeAsVectorFormat( vLayer,
        ShapePath + vLayer.name() + ".shp", "utf-8",
        vLayer.crs(), "ESRI Shapefile" )

    from PyQt4.QtGui import QMessageBox
    QMessageBox.information(None, "MapInfo Mif/Mid Exported to", "" + unicode(MifPath))
    QMessageBox.information(None, "ESRI Shapefile to", "" + unicode(ShapePath))


texts: Defining projection of raster file which lacks prj file

texts:
I have two raster data. They must cover each other. My problem is, one of them doesen't have prj file, therefore I don't know what is the projection of this one. I can not cover this raster with the other one. Of course, I tried define projection function several times, but I don't know which projection should I use.

I use ArcMap 10.6.1


texts: Determining matrix distances between points from different layers along line using QGIS?

texts:
I have been looking for how to obtain a matrix of distances between points through a network with QGIS 3.0.3.

I have only found how to obtain a matrix of linear distances, that option is easy to find wiht the basic QGIS options. However, I can not find the way to extract the example I need: to produce a matrix of distances between points of one layer and others of another, through a linear vector, specifically, the distance of several sampling points to several reservoirs through the fluvial network.

I found it in ArcGIS, but I no longer have that program,

Would anyone know it in QGIS?


texts: Can GDAL fill elevation below a certain value?

texts:
I have several GeoTIFFs that I'm working with to create contours using gdal_contour. One of them crashes the utility. gdalinfo tells me that the computed min values for band 1 and 2 are -340282346638528860000000000000000000000...  This seems to be creating the problem as none of the others are so irrational (think in the single digit range / near sea-level).

Is there a way to "fix" this geotiff by replacing anything below 0 with a reasonable value such as 1 or even 0? And can I simultaneously replace "nodata" values with 0 or 1?


texts: Find 3 closest points in relation to another point

texts:
Using QGIS, I have two shapefiles:

points representing enterprises (about 50k)
points representing towers (about 300)

For each enterprise, I need to figure out the three closest towers and their distances in meters.

Is there a script I can run that can calculate this? Or do I need to use the field calculator/select by location tool?


texts:
I have three buffer radius sizes

radius<-c(1000,700,500)

that I have placed around an array of points (two example points are given here in the file SourcePlace).

What I'd like to do is visualise the three buffer sizes, in different colours, on the same plot. I attach the file with the point information, then

install.packages("sp")
library(sp)
PlaceCoord <- SpatialPointsDataFrame(data.frame(UTM.Easting,UTM.Northing),data=data.frame(PlaceID,Name,lngYear,IDgroup),proj4string=CRS("+init=epsg:32717"))
ptbufN <- list()
for (i in 1:3) {
  ptbufN[i] <- gBuffer(PlaceCoord, width = radius[i], quadsegs = 10)
}
plot(ptbufN[[1]], density=NA, col=1)
plot(ptbufN[[2]], density=NA, col=2)
plot(ptbufN[[3]], density=NA, col=3)
points(PlaceCoord, pch = 20, col = 'red')

What I can't seem to fix is that my polygons are drawn on new plots each time. Maybe there's code to add polygons to existing plots somewhere in ggplot2?