texts: How to split several lines by one polygon and get the distance of each part using sf?

texts:
The question "How to split one line by several polygons" already has an answer here. My question is similar but different: How to split several lines by one polygon?
The context: I have a border and routes followed by vehicles. I want to

decompose the routes, cutting them everytime they cross a border, and
for each part, know what are the distances performed.

Some routes have overlapping segments: the go somewhere and come back using the same segment.
I would like to stay with sf for this problem, to limit the number of packages I'm using in my project.
** Input data**: Here is a simple example with a square as the border and five routes
library(sf)

# Defining the "border", a polygon
four_corners <- matrix(c(0, 0, 4, 0, 4, 4, 0, 4, 0, 0),
                       ncol=2, byrow=TRUE)
pts <- list(four_corners)
my_polygon <- st_polygon(pts)
my_polygon_as_a_linestring <- st_cast(my_polygon, "LINESTRING", group_or_split=FALSE)

# Defining the first "route", with 1 intersection
two_points_of_linestring1 <- matrix(c(0, -1, 2, 1),
                                    ncol=2, byrow=TRUE)
my_linestring1 <- st_linestring(two_points_of_linestring1)

# Defining the second "route", with 2 intersections
two_points_of_linestring2 <- matrix(c(1, -1, 5, 3),
                                    ncol=2, byrow=TRUE)
my_linestring2 <- st_linestring(two_points_of_linestring2)

# Defining the third "route", internal but with overlapping segments
five_internal_points <- matrix(c(1, 2, 2, 2, 2, 3, 2, 2, 3, 2),
                      ncol=2, byrow=TRUE)
my_linestring3 <- st_linestring(five_internal_points)

# Defining the fourth "route", with 3 intersections
four_points <- matrix(c(0.5, 0.5, -0.5, 1.5, 0.5, 2.5, -0.5, 3.5),
                      ncol=2, byrow=TRUE)
my_linestring4 <- st_linestring(four_points)

# Defining the fifth "route", with 5 intersections
five_points <- matrix(c(0.25, 3.5, 1.25, 4.5, 2.25, 3.5, 3.25, 4.5, 4.25, 3.5),
                      ncol=2, byrow=TRUE)
my_linestring5 <- st_linestring(five_points)

# Group the five routes in a spatial dataframe with ID
g <- st_sfc(my_linestring1, my_linestring2, my_linestring3, my_linestring4, my_linestring5)
df_routes <- data.frame(id_column = 1:5)
st_geometry(df_routes) <- g

The geodataframe is thus:
  id_column                       geometry
1         1         LINESTRING (0 -1, 2 1)
2         2         LINESTRING (1 -1, 5 3)
3         3 LINESTRING (1 2, 2 2, 2 3, ...
4         4 LINESTRING (0.5 0.5, -0.5 1...
5         5 LINESTRING (0.25 3.5, 1.25 ...

The expected output is:
  id_column distPart1 distPart2 distPart3 distPart4 distPart5 distPart6 distPart7 distPart8
1         1    1.4...    1.4...        NA        NA        NA        NA        NA        NA
2         2    1.4...    2.8...    1.4...        NA        NA        NA        NA        NA
3         3         4        NA        NA        NA        NA        NA        NA        NA
4         4    0.7...    0.7...    0.7...    0.7...    0.7...    0.7...        NA        NA
5         5    0.7...    0.7...    0.7...    0.7...    0.7...    0.7...    0.7...    0.7...

By using this approach, we decompose all routes by all intersections, which is not correct.
Using st_split does not work with overlapping segments, as described here.
How can I decompose each route by its own intersections?


texts: Error when displaying .shp from postgis with mapserver

texts:
I have postgis database with lines in .shp format in 3D. I made connection in map file:

     LAYER
 NAME paths_3D
 CONNECTIONTYPE POSTGIS
 CONNECTION "host=www.*****.com port=5433 dbname=mountain user=**** password=****"
 DATA "the_geom from paths_3D"
 STATUS ON
 TYPE LINE
 CLASS
    NAME "all_paths"
    STYLE
       COLOR 188 143 143
    END
END
END

Then, I added checkbox in html file pointing on that .shp in postgis database:

<input type="checkbox" onClick="UpdateMap()" name="layer" value="paths_3D" [paths_3D_check] > Paths<br/>

So, when I open my web page and click on checkbox named Paths, it gives me this error:

writeLog(): Unable to access file. mslog.log msDrawMap(): Image handling error. Failed to draw layer named 'paths_3D'. msPostGISLayerWhichShapes(): Query error. Error (ERROR: column "oid" does not exist LINE 1: ...ary(ST_Force_2D("the_geom"),'NDR'),'hex') as geom,"oid" from... ^ ) executing query: select encode(ST_AsBinary(ST_Force_2D("the_geom"),'NDR'),'hex') as geom,"oid" from paths_3D where the_geom && GeomFromText('POLYGON((561889.240401714 5072468.31385521,561889.240401714 5093238.31385536,589594.131887715 5093238.31385536,589594.131887715 5072468.31385521,561889.240401714 5072468.31385521))',find_srid('','paths_3D','the_geom'))

I don't know what is wrong, please help.


texts: Seeking Mobile App to insert customized roads for navigation?

texts:
I want to know if there is a mobile app (for cellphone) where I can insert my line features (roads) and polygons features (fields), all in vector form (.shp) and somehow navigate with that information.
Instead of using downloaded OpenStreetMap maps I want to download my own shapes.

The purpose of that is inserting private roads and to navigate from Point A to point B just as a regular Navigator would do.
The idea is knowing what the fastest route is, the shortest route, the better for trucks, cars, bikes, walking....

The current option is using Navigator from ESRI, but every license is really expensive, and I am looking for a free App.


texts: Merging shapefiles with DBF file, error with duplicate values

texts:
When I try to merge my household data file in a .dbf format with my shapefile by using key variable, I get the following error message:

  Chosen table merge key field IL contains undefined or duplicate values. Key fields must contain valid unique values.

Because, in my dbf file, provinces have several observations. However, my shapefile has just name of the provinces. So how can I merge these two files in GeoDa?


texts: QGIS 2.0, SpatiaLite: labels disappear after setting data defined position (x,y)

texts:
I added to my SpatiaLite table (.osm-import with the "standard" QGIS-OSM-Plugin) two empty columns x and y (double) to define manually the position of my labels (data defined labeling). But after telling the layer about the x and y-column, all labels just disapears!

I tested it in different CRS and with different data types and lengths. I created the columns in QGIS (When I create them in DB-Manager or Qspatialite, they wont appear in the QGIS attribute table)

Im absolutely new to SQL databases!

Thats my table:

UPDATE:

Creating a column with

ALTER TABLE "pelepo_300114_place" ADD COLUMN label_x integer

works but is NOT displayed in Qgis. (I tryed it with integer, numeric, float...)


texts: Data Driven Page Extent based on Polygon

texts: I would like to let the polygons in my Index feature class set the extent of each page,  I am sure this is possible, and the answer is simple, but I can't seem to find the instructions for doing this.


texts: Exracting data fromSpatialPolygonDataFrame conditionally in R?

texts:
For some research I am working on, I have a dataframe called charters of charter schools in New Jersey. Simplified it looks like this:

      school.name          school.id        district.id     lat     lon
oceanside charter school     008229           3400011     39.3635  -74.4350
discovery charter school     228881           3400020     40.7343  -74.1745

I also have a SpatialPolygonsDataFrame of school districts in New Jersey which I got through the Tigris package: nj.school.districts <- tigris::school_districts("NJ"). For regular, non-charter schools, their district.id matches their geographic district.id, as one might expect. However, for charter schools a unique district.id is assigned that doesn't reflect their geographic location at all.

My goal is to assign a geographic.district.id to each school in charters that reflects the district it is actually located in. The nj.school.districts spdf looks like this:

So there are 342 districts, each with a unique GEOID. It is this GEOID that should be assigned to each charter school that falls within the district. I've mapped the charters over the school districts so could do it manually that way, but in reality it would be too time consuming.

I've been able to extract the coordinates of a particular polygon, and to confirm whether or not any schools are in it through this code:

nj.polygon1.coords <- as.data.frame(nj.school.districts@polygons[[1]]@Polygons[[1]]@coords)
%>% select(lat = V2, lon = V1)

xp = as.vector(nj.polygon1.coords$lon)
yp = as.vector(nj.polygon1.coords$lat)

in.poly <- pracma::inpolygon(charters$lon, charters$lat, xp, yp)

but I don't know how to do it at scale or how to then assign the new GEOID accordingly.


texts:
I have mapbasic code that manipulates the obj-column in a point-table. Some rows doesn't have any geometry so the obj-column is empty.  My question is how to check if  there is no object?

All the methods i found in the mapbasic-reference that givs information about the obj-column, like ObjectInfo(), try to access the object which givs an error since there is none.

Is there any way to use SQL to select the rows with non empty obj-columns?

Regards


texts: Loading pre-labeled *.adf raster coverage into QGIS

texts:
I am attempting to load this dataset into QGIS. Following Loading *.adf files into QGIS, I can load the data by opening the hdr.adf file in the Data Source Manager. However, it imports without any labels on the discrete values in the color ramp:

I know the labels do reside in the vat.adf file.
How do I load those labels into the QGIS color ramp?


texts:
I created a line by importing a delimited CSV file and the line appears fine. I can add a label on the line by using the "begin" or "end" attributes but not the "Name" attribute. It is not showing a label. So I assume the value of "Name" is somehow blank although I have it defined in the CSV file:

I define the value as "Name" on the Layer Properties/Label tab:

Is there a way of looking at the attribute values after the import?
After fixing the CSV import file, I get:

I think Stu Smith is right that I was labeling the points and not the line. But how do I do that when the there is no attribute on the line style:


texts:
I have a shapefile with polygons. And I have a global raster file. I want to overlay the shapefile's polygons onto the raster grid and calculate the mean raster value for each polygon.

How can I do this using GDAL, writing the results to the shapefile?


texts:
I am trying to merge two polygon shapefiles in QGIS 2.6.1.

One is older one with all state boundary, which I have (.shp, .shx, .prj, .dbf).
New one is downloaded from WEOGEO (.shp, .shx, .prj, .dbf )
and Coordinate ref system EPSG:4326 - WGS 84 - both.

In the newly merged shapefile the Polygon overlaps. If I set the transparency to high layer below is visible.

I thought of editing attribute table but, both have different attribute fields.

The old one have field like ID, NAME, TYPE, LENGTH, AREA.
But the new one OSM_ID, BOUNDARY, NAME, ADMIN_LEVEL

I am new to QGIS, how can I get the area and length of new polygon if I have to match attribute field of old.


texts:
I have a Lidar Scan from the inside of a building, in the program Leica Cyclon (.ptg .ptz .xyz) and/or Autodesk ReCap (.rcp .pts .rcs .e57). I need to import this data into ArcGIS Scene for 3D viewing.

Is the best way to go about this to use a 3rd party to convert the data into a .LAS and then bring it into ArcGIS Pro?

I used the ASCII 3D to Feature Class tool and then used the Raster Interpolation tool to interpolate multipoints to a raster but it would not work.


texts: Executing function on start up/ on started ArcMap in ArcObjects?

texts:
I need execute a function that init a timer when ArcMap application start.
My extension class:

namespace GMREAddin.Modules.Toolbar
{
    public sealed class eManager : ESRI.ArcGIS.Desktop.AddIns.Extension
    {
       protected override void OnStartup()
       {
            MessageBox.Show("opened from onstartup.");
            initTimer();
       }
    }
}

And have this attribute in the configuration

When I run the application, it doesn't display the message, and does not start the message.

I am using version 10.4.


texts: Print composer will not show maps when reopened

texts:
I completed a print composer project, locked layer and layer style files and locked the map, then saved the print composer, then saved the complete project.  I closed QGIS and when I opened the project again, the saved print composer didn't display the map. All other items, legend labels etc are saved and displayed.

I tried to refresh render/ rectangle/ cache and it didn't work. I had to make the map again.


texts: Geotransformation for polar stereographic?

texts:
I'm currently working to import CANGRID climate data (provided as Surfer Grid ascii, ".grd" files) into ArcGIS. The grid is 95 rows by 125 columns in size.Metadata provides lat/lon of origin (lower left corner), cell size (50km) as well as notes projection as polar stereographic with central meridian (110 degrees W) and latitude of origin (60 degrees N).

After first attempting to convert the .grd to rasters as .ascii and .flt unsuccessfully, I've managed to use GDAL to set the extent and projection, however the dataset does not correctly align with the boundaries of the intended area. See below image.

Is there an accepted geotransformation for polar stereographic that could explain this lack of alignment?

For instance, is there a specific conversion factor or rotation that I should be using?

An example file from the dataset is here: "t201113.grd"

Here's the code I am currently using in GDAL

ds = gdal.Open("t201113.grd")
array = ds.ReadAsArray()

x_rotation = 0
y_rotation = 0
xres = 1
yres = -1

llx = -129.8530
lly = 40.0451
ulx = -175.144
uly = 71.385

input_osr = osr.SpatialReference()
input_osr.ImportFromWkt(ds.GetProjection())

wgs84_osr = osr.SpatialReference()
wgs84_osr.ImportFromEPSG(4326)

wgs_to_nps_trans = osr.CoordinateTransformation(wgs84_osr, input_osr)
x, y, z = wgs_to_nps_trans.TransformPoint(llx,lly)

geo_transform = [ x, xres, x_rotation, y, y_rotation, yres ]

ncol = ds.RasterXSize
nrow = ds.RasterYSize

out_driver = gdal.GetDriverByName("HFA")
out_ds = out_driver.Create(t201113.img", ncol, nrow, 1, gdal.GDT_Float32)

out_ds.SetGeoTransform(geo_transform)

out_prj = 'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-110.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",60.0],UNIT["50_Kilometers",50000.0]]'

out_ds.SetProjection(out_prj)

out_ds.GetRasterBand(1).WriteArray(array)
out_ds.GetRasterBand(1).SetNoDataValue(1.70141e+038)
out_ds.FlushCache()
out_ds = None
`

Also, here's the projection info, as defined by the input, i.e. from "GetProjection()":

  'PROJCS["North_Pole_Stereographic",GEOGCS["GCS_WGS_1984",DATUM["WGS_1984",SPHEROID["WGS_1984",6378137.0,298.257223563]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Stereographic"],PARAMETER["False_Easting",0.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",0.0],PARAMETER["Scale_Factor",1.0],PARAMETER["Latitude_Of_Origin",90.0],UNIT["50_Kilometers",50000.0]]'

And the input GeoTransform:

(-0.5, 1.0, 0.0, 94.5, 0.0, -1.0)

Lat, longs of the grid coordinates are also provided, and when view in the projected coordinate system look like below. When the geotransform is defined by coordinates of the lower left (yellow) or upper right (pink) cordinate, I can effectively set the extent, but there remains alignment issues throughout the raster.


texts: How to have a Raster output of a ppm model in spastat

texts:
I am trying to compare the outputs of species distribution models. My current code uses the following function to calculate TSS scores.
TSS_Score<-function(TrueBinary, Raster) {

  Raster[]<-as.numeric(Raster[] > 0.5)

  ConfusionMatrix<-crosstab(TrueBinary,Raster)

  a<-ConfusionMatrix[2,2]

  b<-ConfusionMatrix[2,1]

  c<-ConfusionMatrix[1,2]

  d<-ConfusionMatrix[1,1]

  Sensitivity<-a/(a+c)

  Specificity<-d/(b+d)

  TSS<-Sensitivity+Specificity-1

  return(TSS)
}

This code works but I need a raster as the model output for it to work. I have run a ppm model from spatstat and then used the predict.ppm function. However the output of predict.ppm is an image class.
Is there anyway to get a Raster as a predict output?
Or to convert the image into a raster?
EDIT: Here is the code for my ppm model. The covariates are images files as the function doesn't seem to take rasters?

PPMmodel<-ppm(PointsPPP, ~slope, covariates = list(slope =
AMtemp,MDRtemp,Isotherm,AMprecip,SeasonPrecip, TreeCover50))
predictPPM<-predict(PPMmodel)
predictPPM
real-valued pixel image
128 x 128 pixel array (ny, nx)
enclosing rectangle: [92.192, 141.01] x [-10.999, 28.549] units


texts: Select specific attribute from SDE.ST_Geometry UDT? (without using accessor functions)

texts:
Oracle 18c; ArcGIS 10.7.1

I have an SDE.ST_Geometry polyline feature class.
As an experiment/learning exercise, I am trying to learn how to interact directly with the user-defined type (UDT) SHAPE column.
I can select the SHAPE column in an SQL client which returns a collection of attributes from the UDT:
select
    shape
from
    a_test_fc

SHAPE
--------------------------------------------------------------------------------
ST_GEOMETRY(4,37,671828.7852,4861707.1775,671839.2487,4861717.6411,NULL,NULL,NULL,NULL,0,32.9250330071201,26917,(BLOB))

1 rows selected.

And I can see that the docs describe the UDT attributes as follows:

The attributes of the spatial type represent the following
information:

Entity: The type of geometric feature stored in the spatial column (linestring, multilinestring, multipoint, multipolygon, point, or
polygon), the value of which is a bit mask derived from the
st_geom_util stored procedure.
Numpts: The number of points defining the geometry; for multipart geometries, this includes the separators between each part, one point
for each separator.
Minx, miny, maxx, maxy: The spatial envelope of the geometry
Area: The area of the geometry
Len: The perimeter length of the geometry
SRID: Contains the identifier for the geometry that links it to its associated spatial reference (coordinate system) record in the
ST_Spatial_References table
Points: Contains the byte stream of the point coordinates that define the geometry

Question:
Is there a way to select specific attributes from the ST_Geomtery UDT?
I'm aware that there are "accessor functions that take a single ST_Geometry as input and return the requested property value as a number." So there isn't a practical need to directly select individual attributes from the UDT — I could just use those accessor functions if required. I'd just like to determine if it's possible to interact directly with the UDT attributes as a way to learn more about spatial database development.
I tried researching non-GIS UDTs, and how to select their attributes (example: I tried selecting SHAPE.POINTS). But that just resulted in an error: ORA-00904: "SHAPE"."POINTS": invalid identifier. And yes, I'm aware that the POINTS attribute is stored as a byte stream and won't be useful in any practical sense.


texts: How can a gui be visible till all the work is done? (QGIS plugin, python)

texts:
I have made a gui in QTDesigner for QGIS2.10 (as a dialog).
When I press the ok button (if result == QtGui.QFileDialog.Accepted:) everything works just fine but the gui is disappeared.
How can the gui stay there - be visible till all the work is done (it is set to be on top)?
Do I need something like a timer or should I look for some kind of relative property?


texts: Is there a way to calculate length of linestring from one point to another?

texts:
If I have one linestring and coordinates of two points on that linestring, what is the easiest way to calculate not shortest distance from those two points, but distance over the given road  (given linestring) in PostGIS?

So in the given example given linestrin is whole line and I want to calculate length of red line between points A and B.
Any help with this?


texts: How to display hidden line types properly when printing from AutoCAD?

texts: I want to show the pline type in hidden line type but when I change it I can't show it in printed paper, if I explode it the line are so small that can't show at print, so I want to change the pline to single line, what should I do?


texts: SelectLayerByAttribute_management searching with different datatypes

texts:
I am trying to use the SelectLayerByAttribute_management in a tool.
The tool works like that:
The User puts in the Feature which Layer should be selected by attribute -> "inputFeature".
Then he chooses the Field of the table
-> "inputFieldname".
At the end he types in the value which will highlight the Layers by Attribute -> "inputValue".
That is my code, which works just with numeric value.

def execute(self, parameters, messages):

    """The source code of the tool."""

    inputFeature = parameters[0].valueAsText
    ap.AddMessage("Feature: "+ inputFeature)
    inputFieldname = parameters[1].valueAsText
    ap.AddMessage("Fieldname: "+ inputFieldname)
    inputValue = parameters[2].valueAsText
    ap.AddMessage("Value: "+ inputValue)

    whereClause = """{} = '{}'""".format(ap.AddFieldDelimiters(inputFeature, inputFieldname), inputValue)

    ap.SelectLayerByAttribute_management(inputFeature,"NEW_SELECTION", whereClause)
    return

If I choose a field where just Strings are stored and I type in for the inputValue a String, I get this Error:

  ExecuteError: ERROR 000358: Invalid expression. Failed to execute
  (SelectLayerByAttribute).

If I choose a field where just Numbers are stored (doesn't matter which type...double, float, integer) and I type in for the inputValue a Number, I get results. Or better said, all values which have been found are highlighted with SelectLayerByAttribute.

So what is the solution for typing in any datatype I want...?


texts: OSGEO4W Checking GDAL version with `gdalinfo --version` returns nothing

texts:
I am using OSGEO4W 32-bit on Windows 10.  I am trying to check the version of GDAL.
 Following the instructions here, when I input gdalinfo --version in the OSGEO4W shell, I get just a blank line.

Why am I receiving just a blank line and not the version number?


texts: Building QGIS 2.16.3 Plugins on Mac OS X Sierra

texts:
I'm having difficulty building QGIS plugins using the new Mac OS Sierra.  I can compile the resources using

pyrcc5 resources.py resources.qrc

and the plugin displays in QGIS, but the message
'This plugin is broken.  No module name PyQt5'

I have PyQt5 installed on my computer, but it seems that it may not be in Qgis 2.16?

I'm also not certain that PyQt4 works in Sierra.  I have some experience with python and building plugins, but am pretty new to command line.


texts: How to print a large map easily from QGIS where I can later zoom in on minute details?

texts:
I have a cadastral map with more than 20000 plots. For a ground survey, we need to use this map to capture details of each plot.

I want to make a pdf file that can opened on mobile phones for the survey. Here is how the map should look like

Actual map on map canvas

This is how the zoomed in version should look like.

This is how much zoomed in level I require. I want it to zoom in in the pdf just like it zooms in on the qgis desktop

But, unfortunately, I need to zoom in on each map tile and cut it, merge in paint and export as pdf. Otherwise, 90% of details are lost by printing it via the print composer. This is an extremely painstaking procedure and I wonder if there is an easier way.

I tried to load it as extremely high res map on print composer and exported as pdf, still the details are missing.

Is there a possibility to print the entire map outside of the map canvas? Then I can zoom the map to the level I want and print the entire map? This process takes me somewhere close to 10 hours and it is precious time lost.


texts: Linking SharePoint database with ArcGIS?

texts:
The organization I am working for has a sharepoint database that multiple people update every time they do an inspection of natural resources on a parcel in the county.

This database has the following fields:

*Parcel Number

*Date

*Notes

*Inspector

*Action

I have created a feature class with a centroid with these attributes placed in each parcel that has been inspected.

Is there a way to dynamically link the sharepoint database to the geodatabase so that when they go in and change/add anything to the notes section it actively changes (and adds new points to newly listed parcels) the attributes in the featureclass?

If not, his layer only needs to be updated monthly. I am able to export it to a excel or csv file. Is there a quick way to create the centroids layer from this table?


texts: Tiger geocoder ask for password for each step while downloading the data

texts:
I'm downloading the data for geocoder with the respective scripts:

nation_script_load.sh
ma_load.sh

The problem I'm having here is that for each of the files I asked to type the password, although the correct password is in the files:

WGETTOOL="/usr/bin/wget"
export PGBIN="/usr/lib/postgresql/9.3/bin"
export PGPORT=5432
export PGHOST=localhost
export PGUSER=postgres
export PGPASSWORD= "my_password"

I have tried by removing the quotes but isn't working either:

export PGPASSWORD= "my_password"
export PGPASSWORD= my_password


texts: Not able to add .jpg image in ArcGIS Online

texts:
I am trying to show an image in the customise popup for every feature in ArcGIS Online. For this purpose I am tring to add .jpg format image to ArcGIS Online. I believe this is an acceptable format but it still not letting me upload the file and is giving me the error that "The file is not supported".

My ArcGIS Online is a public account if it makes a difference and the image I am trying to upload is less than 1MB in size.


texts: Setting Output Extent (Environment Setting) using ArcPy?

texts:
I'm trying to figure out how to set my output extent before creating Thiessen Polygons in a loop of ~100 points shapefile. (ArcMap 10.0 - Python)

for i in range(3,51,1):
    arcpy.CreateThiessenPolygons_analysis("AVF%i" % i,"AVF%iVoronoi" % i,"ONLY_FID")

Basically I want my extent to be the same for each iteration.

I first found:

arcpy.Extent(293490,5898230,316430,5930230)

then someone told me about this ESRI help page:

#use 'CURRENT' if running from ArcMap, when published use MXD on disk

    mxd = arcpy.mapping.MapDocument("CURRENT")

    df = arcpy.mapping.ListDataFrames(mxd)[0]
    newExtent = df.extent
    newExtent.XMin, newExtent.YMin = 293490,5898230
    newExtent.XMax, newExtent.YMax = 316430,5930230
    df.extent = newExtent

but in any case, when I run my loop, the size of the output is not reflecting the changes.

ps: https://stackoverflow.com/questions/8996965/output-extent-environment-setting/8997122#8997122


texts:
I have a point dataset of solar farm locations.  The points are grouped according to which substation they feed into (see below).  I would like to test the hypothesis that electrical output variability within each group is related to the pattern the points form.  That is, 10 points within a 10x10km area may resemble a straight line, wavy line, perimeter of a circle, semi-circle, triangle, dumbbell etc.  However, I’m struggling to find any algorithms which can help to categorise the groups of points or even the correct terminology to describe the possible patterns.

Most spatial analysis techniques simply describe data as “random”, “uniform” or “clustered” but I want to classify the types of cluster.

I’ve found I can identify lines in point patterns with PAST (http://folk.uio.no/ohammer/past).  Hierarchical clustering in R looks suitable for dumbbell i.e. two groups forming a larger group.  But I can’t identify any techniques for the other possible patterns.


texts: How can I automatically split a spatial layer to a grid

texts:
Possible Duplicate:
  Divide vector in equal smaller parts in QGIS or similar

I am a software engineer and the company I work for has a project that needs some spatial data manipulation in order to provide routes to users. The client has provided me the spatial layer through a Microsoft SQL Server 2008, encoded as geometry polygons. The table is as:
zone_id (int), zone_geometry (geometry).

In order to make the routing algorithm more precise, I need to split each zone in tiles, so I thought that I need to divide the each zone to a grid. Is there an automated way to do this?
Are there any references or any wiki for what I need?
I searched the answers in the stack exchange but to be honest I could not understand because I have no previous experience with GIS.

Thank you


texts: Calculating Moran's I for buffers of forest cover?

texts: I'm trying to figure out if buffers I have created in ArcMap (measuring percent forest cover in a 1 km distance from a land use map within the buffer) around points at which I will be collecting presence/absence data of mammals exhibit spatial autocorrelation. Basically I want to ensure each of my buffers of forest cover is a distinct unit of measurement even if they are close together in space. Each buffer measures percent forest cover as the explanatory variable for the response variable of mammal presence/absence gathered at the central point of the buffer. I want to make sure these variables of forest cover in each buffer are actually independent variables that can be used as explanatory variables for each individual buffer (and thereby each individual point where I will measure the response variable, mammal presence/absence) without having spatial autocorrelation. I'm a bit new to this so I'm not sure exactly how this should be done with Moran's I. For example, what the input field for this tool in ArcMap should be? Is there a way I could use the values for percent forest cover for this? I determined the amount of total land use in each buffer with the tabulate area tool, and then actually determined percent forest cover in excel, so I don't know if there's a way to merge the land use values from each buffer derived from the tabulate area tool to the corresponding buffer polygon file so that these can be included into the Moran's I tool?  I'm assuming the input feature class should be all the buffers of interest merged together?


texts: How to convert KML to shapefile in ArcMap 9.3.1

texts: I would like to know about conversion tools available  for KML to shapefile formats. In ArcMap 9.3.1


texts: How can I highlight or color polygons intersecting another layer based on certain values?

texts:
I have one layer of polygons that stretch across the U.S., and another layer of all counties in the U.S. In the counties layer is a column with population density for each county.

What I am trying to do is to filter the data by polygons that intersect with a county that has a specific population density. For example, I want to highlight all polygons that intersect a county that has a pop density of <50. Ideally, it would be nice to have it create a new layer based on these filtered polygons so that I could color them differently.

Is this possible to do in QGIS using an expression of some kind?


texts: Syntax for queryFilter.WhereClause in C#

texts:
I need help with a queryFilter WHERE clause that is giving me a headache.

All it does is to convert the search to uppercase and perform the search whether the result is uppercase or not (e.g. it will find "Jones" or JONES").

I have the following code in my VBA app:

Dim pQF as IQueryFilter
Set pQF = new QueryFilter
pQF.WhereClause = "upper(" & strFieldToSearch &") LIKE '%" & UCase(Replace(strSearch,"'","''")) & "%'"

In my C# code I have tried various ways including:

queryFilter.WhereClause =  searchField + " LIKE '%" + search.Replace("'","''") + "%'";

queryFilter.WhereClause = strSearchField.ToUpper() + " LIKE " + "'%" + strSearch.ToUpper() + "%'";

But it can't find the record. What am I doing wrong?


texts: "Failed to get property: gid" error in geoserver sql view

texts:
I have a postgis data with 777 features. In geoserver, i created a sql view with another table. In the wms request i can see all features, in qgis also i can see all features. but in wfs output as csv or geojson i only get 55 features. On doing wfs with openlayers it gives error

<ows:ExceptionReport version="1.0.0"
  xsi:schemaLocation="http://www.opengis.net/ows http://localhost:8080/geoserver/schemas/ows/1.0.0/owsExceptionReport.xsd"
  xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xmlns:ows="http://www.opengis.net/ows">
  <ows:Exception exceptionCode="NoApplicableCode">
    <ows:ExceptionText>java.lang.RuntimeException: Failed to get property: gid
Failed to get property: gid
java.lang.String cannot be cast to com.vividsolutions.jts.geom.Geometry</ows:ExceptionText>
  </ows:Exception>
</ows:ExceptionReport>


texts:
I want to create an urban center boundary (in ArcGIS v10.3) by using a population-grid dataset.

So, I selected high density cells, that is to say cells with population over 1.500 and hereupon I want to check each of the high density cell if it touches another high density cell, provided that diagonal touching is not allowed.

Any ideas or tools for doing this task quickly?


texts:
I have a geoprocessing service in ArcGIS for Server 10.1.

I would like to share this service with a number of external organisations that also have ArcGIS Server 10.1.
i.e I dont want them hitting my service, but for them to have the ability to easily re-publish the service at their end.

Ideally, I would like to bundle this into a simple installer, which prompts end-user for some admin credentials to ArcGIS Server, and then makes use of a python script to publish the service via the new ArcGIS Server Administration API.

However - To publish a GP service, the help states:

  To publish a geoprocessing service, you need a result in the Results
  window

I therefore can't think how an installer would be able to publish a GP service from a python script, as sounds like it is restricted to being published from within ArcGIS for Desktop.

I also believe that it is not possible to achieve publishing a GP service from a GP package either.

My question is, are my assumptions right, and what would be the best workflow for me to share my GP service, so that my end-users can publish it to a GP service in the simplest approach possible.

i.e. I want to avoid any human error, permission issues, and minimise the time and effort it would take for the end-users to publish the service seamlessly to their ArcGIS for Server.


texts: ArcGIS Javascript API edit data programmatically

texts:
I have a form that auto generates a unique number. I want to update/edit an attribute with this number. Is this possible using one of the Esri editing widgets? Attribute Inspector, Editing widget or Template Picker etc...or would this be easier to implement with an applyEdits() call?

All the user needs to do is drop a point on the map. The only attribute besides the Esri objectId will be this form generated number. Wondering how to tie this auto generated form number to an attribute on the point automatically so that the user does not need to manually edit the attribute.


texts:
I imported my shapefile dat to R with readOGR, coordinates are given in easting and northing, I'd like to change it to latlon and tried:

dat_latlon <- spTransform(dat, CRS("+proj=longlat +datum=WGS84"))

However, I get this error message:

Error in spTransform(xSP, CRSobj, ...) :
  No transformation possible from NA reference system

And str(dat) shows that:

  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slots
  .. .. ..@ projargs: chr NA

So I guess I have to assign the projection first?
Somewhere here (or stackoverflow) I found this:

dat@proj4string <- "+proj=tmerc +lat_0=0 +lon_0=9 +k=1 +x_0=3500000 +y_0=0 +datum=potsdam +units=m +no_defs"

but I only get an error:

Error in checkAtAssignment("SpatialPointsDataFrame", "proj4string", "character") :
  Zuweisung eines Objektes der Klasse “character” ist für @‘proj4string’ in einem Objekt aus Klasse “SpatialPointsDataFrame” nicht zulässig; is(value, "CRS") ist nicht TRUE

(sorry for German - apparently it's not valid to assign an object of class character in an object of class SpatialPointsDataFrame).

Now I am stuck, can you help me?


texts:
I am running ArcSDE 10 with SQL Server 2008 R2 Standard Edition.  I am new to SDE and SQL Server, but I understand that SQL Server has the ability to create relationships between tables and maintain certain referential integrity rules.

ArcGIS has relationship classes that act similarly, but a relationship class does not have all the features of SQL relationships and does not result in a SQL relationship in the ArcSDE database.

Is it possible to create relationship classes in ArcGIS for an ArcSDE database and create relationships for the same tables in SQL Server?  By doing so, I will be able to utilize these relationships whether I am working with the data in ArcGIS or in SQL Server Management Studio.  Will the two types of relationships conflict with each other or otherwise hinder performance?


texts: PostGIS raster not loading in QGIS (version >3.10 on MacOS and Windows)

texts:
When loading PostGIS rasters (through the brower or via DBmanager) from my database into QGIS (LTR versions on both MacOS (Big Sur) and Windows 10), these fail to render. Either nothing shows up, or a band is shown covering most of the UTM zone in which the raster is positioned (see image below. In this case, the original raster covers less than a square kilometer).
In QGIS 3.10 everything works perfectly. In QGIS versions after that, despite the rasters being listed in the layers' panel, nothing is being visualised (or I get the odd error below).
I noticed a couple of posts have mentioned similar issues, but without any response. I'll file an issue in the QGis issue queue, and will update here when I can.


texts: How to load big ASCII rasters in Revolution R?

texts:
Is there anyone in this site using Revolution R with big spatial data?

I am trying to load raster files (ascii format) into revolution R using the xdf format, with no luck. Did anyone tried this already?


texts: Downloading RGB Sentinel-2 with Google Earth Engine and Python

texts:
I am fairly new to the Google Earth Engine platform, and I want to create a dataset for my machine learning project using the aforesaid platform. As far I have this code for downloading a single picture:
import ee

# Trigger the authentication flow.
ee.Authenticate()

# Initialize the library.
ee.Initialize()

geoJSON = ... # coordinates of rectangle
coords = geoJSON['features'][0]['geometry']['coordinates']
aoi = ee.Geometry.Polygon(coords)

ffa_db = ee.Image(ee.ImageCollection('COPERNICUS/S1_GRD')
                       .filterBounds(aoi)
                       .filterDate(ee.Date('2020-08-01'), ee.Date('2020-08-31'))
                       .first()
                       .clip(aoi))

url = ffa_db.select('VV').getThumbURL({'min': -20, 'max': 0})

After that, I just download a file with this command: !wget {url} (From Jupyter Notebook)
The question is fairly simple, but I had much struggled with it due to the lack of Pythonic ways of doing such things. I searched GEE Community tutorials and Github, but have no luck with them.
How can I download RGB images from Sentinel-2 using Python?
Better to download to a local directory rather than a Google disk.


texts: Defining neighborhood settings in Focal Statistics with Calculate Value using ModelBuilder

texts:
I am using ArcGIS Desktop 10.6 with Spatial Analyst.
I use ModelBuilder, and I created a model with focal statistic. There are 2 different workflows now.

calculate the radius value: This value depends on the vertex distance of input- I get this number using with Calculate Value tool
run focal statistic - Neighborhood settings: circle, map unit, where Radius would be equal the calculated value
Now I have to run the first part of the model, than I have to type in the number manually.

Is there any way to use this number automaticly without typing?


texts: QGIS 3.4 - In Layout View can I use the properties of one item as a variable in an expression for another item?

texts:
In my layouts I have a legend and a references text box. I want the references text box to be next to the legend in every layout however the width of the legend changes often. Both items are in the bottom right corner of the layout.

Both items' position reference point is bottom right. To position the references box correctly each time I take the X position of the legend, subtract the width of the legend, and then subtract the distance I want between them. The result is the X position of the references box.

Instead of doing this manually can I use an expression like this for the X position of the references box:

"Legend_X" - "Legend_Width" - 5


texts: pgrouting - topology creation with explicit node file

texts:
Learning pgrouting, I would like to build a simple public transport routing engine.

Starting with a simple metro network (DC in my case, see excerpt below), I would like to combine the line geometries (LineStrings) and station geometries (Points) to produce a routable network following projects like these from @underdark!

What is not clear to me is if there is a way to utilize explicit station geometries as input to the pgr_createTopology function to create the network topolgy?

A screenshot of lines geometry overlayed with station locations: .


texts: Finding nearest point outside series of geometries using PostGIS?

texts:
I have a PostGIS database with a very large amount of geometries (circles), most of them overlapping.
I'm trying to come up with a way to find the nearest point (from your current location) that is not within any of these circles.

Is this possible with a direct query?

So far I thought of a fishnet approach of 1x1 km cells that check for overlap with any of these circles. This works, but is not very accurate.


texts:
Using Arc 10.3.1 ModelBuilder, I need to populate a .dbf field (pfo_pct) with a variable value using the Calculate Field tool. However, the .dbf is an empty table (no records) and I get a Warning 000405: No records within table.  The same problem happens outside of ModelBuilder when I try to Field Calculate on an empty table.

I believe that the problem could be resolved if I could insert an empty record into the .dbf prior to the field calculation, but I don't see any tools that accomplish that, nor any environment settings that would help.

Sadly, I am not proficient with Python...


texts: Seeking PLSS Quarter of Quarter Sections in Oregon

texts:
I'm looking for a quarter of a quarter of a section (PLSS) which is 40 acres in size for the state of Oregon. I looked at several government and Oregon state web site and could not find it. The only thing I can find is Township, Range, Section.

Earlier threat (2 years ago) refers to creating a quarter of a quarter grid using QGIS, but I'm using ArcMap 10.5 and looking for a completed shapefile, not creating one from scratch.

See image bellow (circled in red).


texts: arcmap export layer as a image with different scale

texts:
In our web application, we load too many layers cause the map loading very slow. I am thinking  of export some static layers as images with different scale. There are totally 10 scales in our application, hence I want to export 10 images for a static layer (each scale 1 image). When the application load map data, I only load the image file for the static layer. When the map zoom, I will change different image according to the zooming scale.

I am using arcmap 10.2.

I feel it can be performed, but I haven't get the way to do it. Please advise me, thanks.


texts: Canonical form for romanizations of Arabic town names in Syria mapping to latitude/longitude

texts:
Is there a Python package or database source which would allow me to go from various romanizations of Arabic town names in Syria and Iraq to latitude/longitude pair?

Here is an example where the left romanization is one kind used by a US Government record collection and the other appears in a Syrian government press release.  The Government names seem to be more well known and look-upable than the names used by the Syrian Government.  I would like to find a canonical reduction of the two into a single reliable name:

Al Bāqāt::al-Bakat
Al Ismā‘ilīyah::al-Ismailiya
Arjul::Arjal
Batrānah::Batrana
Bīshah::Bisha
Buraydah::Barda
Burj ‘Izzāwī::Burj Ghazawi
Burj Subaynah::Burj Sabna
Buslah::Bassila
Buwayḑah al Kabīrah::Bouaida Kabira
Buwayḑah aş Şaghīrah::Bouaida saghira
Ghurayrīfah::Gharirfa
Jabal al Mudawwar::Jabal al-Madwar
Jafr Manşūr::Ja’far Mansour
Jubb al Khāfī::Jub al-Khafi
Jubb at Tīnah::Jub al-Tina
Judaydah::Jadida
Kafr Ḩawt::Kafro Hoot
Khanāşir::Khanaser
Khirbat al Ma‘ājīr::Khirbet al-Muaijer
Madāyin Kabīrah::Madaen Kabir
Munbaţiḩ::Munbateh
Rasm al Bās::Rasm al-Bassas
Rasm al Karkūr::Rasm al-Karkour
Rasm Shawkān::Rasm Shoukan
Rujaylah::Rjaila
Şuwayyān::Soyan


texts: Raster Properties to TextFile using ArcPy?

texts:
I have multiple rasters that I want to get the mean of the all the cells. This is what I have:

env.workspace = r"C:\Users\pyb609\Desktop\Updated
Files\GIS\Dark\Shapefiles\15M"
outfile = open(r"C:\Users\pyb609\Desktop\UpdatedFiles\GIS\Dark\Shapefiles\15M\15M_Mean.txt", 'w')

#Check out the Arcgis Spatial Anaylst extension licence
arcpy.CheckOutExtension("Spatial")

# Get a list of Rasters
rasterList = arcpy.ListDatasets("*","Raster")

for inRaster in rasterList:
    # just to view progress
    print inRaster
    outfile.write(inRaster)

    Mean =arcpy.GetRasterProperties_management(inRaster, "MEAN")
    textinfo =(inRaster, Mean)
    textstring =str(textinfo)
    outfile.write(textstring)
    outfile.close()

It gives me the first file but then I get an error:

Traceback (most recent call last):
File "C:\Users\pyb609\Desktop\Updated Files\GIS\Dark\MeanRaster.py", line 26, in <module> outfile.write(inRaster)

ValueError: I/O operation on closed file


texts: Hotlink vector/raster files to open in current workspace in ArcGIS 10 similar to MapInfo

texts:
I want to add datasets to ArcMap and place them in a particular group layer by simply clicking on hotlinked polygons in a shapefile. The same process is easy in MapInfo and I want to implement the same for ArcGIS 10.

I found the following post http://forums.arcgis.com/threads/8236-Hyperlinks-in-ArcGIS-10 that describes the process in 9.3.1

On properties-->html popup you can load url and any file will open in an external program but I want it to add to the current workspace.

A solution similar to mapinfo--layer control--hotlinks will be ideal but a python based solution will be acceptable as well if it can be added to the relevant vector layers within an mxd.


texts: Providing a qualitative description of the results from DBSCAN clustering

texts:
I have run a DBSCAN clustering algorithm in QGIS 3 on bus stops in my specific study area where my minimum cluster size is 50 and maximum distance is circa 9km.

After this I created a concave hull of the clusters and found the centroid of each of these.

How can I qualitatively define that centroid in the context of this methodology that is easy to understand for a non-technical reader?


texts: Writing "m" after contour line label numbers

texts:
does anyone know what is the expression to display the letter "m" after contour line numbers in Arcgis?

I want it to be like this.

I believe is in this window wizard, just don't know what expression to use


texts: Seeking world map with boundaries of regions and communities for commercial use?

texts:
I am looking for a vector file worldmap which should include as many boundaries reaching to a very small level (like ZIP-Codes/small communities/etc) as possible.

The data of Natural Earth is quite good, but does not go that much in detail. The GADM data would be nearly perfect, but it is not allowed for company use.

Requirements:

Consistent boundaries for nearby countries (no overlapping!).
As much detail as possible (refering to boundaries).
Country boundaries and subdivisions should be separated.
The data is for company use, therefore I should be able to buy it, or it should be free to use for companies.

I also thought of extracting it from the OpenStreetMap data. Do you think this is possible?


texts:
Does ArcGIS Pro afford the ability to batch create chart elements?
I am aware of how to create a single chart by selecting the table in the TOC, right click and "create chart". I can select the "series" group that I want to show in that chart (such that I can hide the lines for series groups that I do not want shown).
My table has dozens of sites with multiple test results and test dates.
My need is to produce the same format chart for each of my sites. Only showing data for a single site on each chart. I can then manually add the chart(s) to the layout as needed.
I'm picturing a GP tool to create chart element.


texts: Trying to use QGIS to make AK election results maps

texts:
I have QGIS.

I have a current shapefile of the political boundaries of Alaska's boroughs and census areas, acquired from the state's GIS engineers.

I have a PostgreSQL DB of AK election results going back to 1996 that I have beaten until it produces results by borough and census area.

All I want to do is make pretty little results maps that paint the political subdivisions degrees of red and blue.

I have absolutely no idea where to start or where to look for guidance.  Can anyone point me in the direction of a "QGIS for Drooling Imbeciles" type of guide?

I had previously been using the MapViewer 7 and MapViewer 8 software, but this poses problems because the AK gsm sample files that come with the software are relentlessly inaccurate.


texts: How can I programmatically open an OnPoint map to specific coordinates?

texts:
We have an OnPoint server at work and I'm wondering how to open the mapping service and having it zoom in to my specified latitude and longitude.

The reason I ask is because we have an application at work that requires our users to manually type in the latitude and longitude into our system for searching. Our application does have this information so we'd like to offer the users a link.

Here is an example URL that you could plug into Google Maps.

Here is an example of a GIS system running OnPoint.

I do see if their promotional PDF clearly states that you can perform searches via URL queries but I see zero documentation.

 Launch OnPoint from external applications by sending URL to OnPoint Website


texts: Determining lat/long of location when in UTM projected map?

texts:
I am working in an ArcMap MXD, where all my layers and the dataframe are in UTM 19N.  What's the easiest method to determine the lat/long of any point (a random point of interest) within my map?

Do I have to re-project a layer to do so?  Or is there a method where I can just click and get a value in lat/long?


texts: Is it possible to "crack open" the QGIS GUI to see what coding commands are executed whilst working?

texts:
I would like to see the coding that is executed when working in QGIS e.g. using the tool bars and performing tasks in the GUI. Is there a way to see what coding is going on behind the GUI?

I would like to slightly edit what some of the functions do and need to see the code executed when I use the functions.


texts: How to render list of static images with coordinates in OpenLayers

texts:
I am novice in OpenLayers. Have got list of images and their coordinates (in another list) extracted from  big TIFF file. Requirement is to rendered the images with coordinates to get same image as in TIFF file.
Have tried the display of static image example present in OpenLayer site and able to rendered it with single image. When I tried similarly with multiple images it is getting rendered on top of previous one.
var projection = new Projection({
  code: 'xkcd-image',
  units: 'pixels',
  extent: extent,
});

var map = new Map({
  layers: [
    new ImageLayer({
      source: new Static({
        url: images[0],
        projection: projection,
        imageExtent: extent,
      }),
      opacity: 1,
      zIndex: 1,
    }),
    new ImageLayer({
      source: new Static({
        url: images[3],
        projection: projection,
        imageExtent: extent,
      }),
      opacity: 1,
      zIndex: 4,
    }),
    new ImageLayer({
      source: new Static({
        url: images[4],
        projection: projection,
        imageExtent: extent,
      }),
      opacity: 1,
      zIndex: 100,
    }),

    new ImageLayer({
      source: new Static({
        url: images[6],
        projection: projection,
        imageExtent: extent,
      }),
      opacity: 1,
      zIndex: 100,
    }),
    new ImageLayer({
      source: new Static({
        url: images[7],
        projection: projection,
        imageExtent: extent,

      }),
      opacity: 1,
    })
  ],
  target: 'map',
  view: new View({
    projection: projection,
    center: getCenter(extent),
    zoom: 1,
    maxZoom: 8,
  }),
});

I am not able to find any clue how to rendered the images based on their coordinates to get the image correctly. Any suggestion (with example)? Using OpenLayers v6.


texts: GeoServer CSS rotation with env and property name

texts:
Using GeoServer CSS extension

 mark-rotation:  [env('rot',2)];

ends in

    <sld:Rotation>
        <ogc:Function name="env">
            <ogc:Literal>rot</ogc:Literal>
            <ogc:Literal>2</ogc:Literal>
         </ogc:Function>
     </sld:Rotation>

How should my style look like to add a property name to rotation:

 <sld:Rotation>
    <ogc:Add>
    <ogc:PropertyName>mystyle</ogc:PropertyName>
       <ogc:Function name="env">
          <ogc:Literal>rot</ogc:Literal>
          <ogc:Literal>2</ogc:Literal>
       </ogc:Function>
    </ogc:Add>
  </sld:Rotation>

I have tried

 mark-rotation:  [env('rot',2),mystyle];
 mark-rotation:  [mystyle],[env('rot',2)];


texts: Calculate MSAVI (Modified Soil-adjusted Vegetation Index) in Google Earth Engine

texts:
How do you write the MSAVI formula in JavaScript function for Google Earth Engine?

I am trying to calculate the MSAVI as a function. The aim is to map the result over a Sentinel-2 composite and not a single image.

I tried the same thing for the OSAVI Index and it worked fine:
OSAVI = (NIR-Red)/ (NIR+Red+0.16)

function addOSAVI(input) {
  var nir = input.select(['B8']);
  var red = input.select(['B4']);
  var osavi = nir.subtract(red).divide(nir.add(red).add(0.16)).rename('OSAVI');
return input.addBands(osavi);

}

I need a function in the same manner as the OSAVI function, but I have problem with implementing sqrt() and pow() correctly.

function addMSAVI(input) {
  var nir = input.select(['B8']);
  var red = input.select(['B4']);
  var msavi = (nir.add(1).multiply(2).subtract(sqrt(pow(nir.add(1).multiply(2),2)).subtract((nir.subtract(red)).multiply(2)))).divide(2).rename('MSAVI');
return input.addBands(msavi);

}

More information about MSAVI:
http://wiki.landscapetoolbox.org/doku.php/remote_sensing_methods:modified_soil-adjusted_vegetation_index


texts: Overpass API : add center to custom output

texts:
I'm trying to get a list of neighbourhoods with the latitude and longitude + the admin div where each is located. It works fine with nodes, however I would like to do the same with relations ( in that case I just want the center ) but I can't find a way to do this.
I know about "out center", but I don't understand how I can associate it with the "make out" statement, I imagine it's possible to create the centers of relation this way and then add to my final output?
This is what I have so far, probably not efficient at all
    [out:csv(::lat,::lon,_row;false)]
    [timeout:18000];
    {{geocodeArea:"Estonia"}}->.searchArea;
    (
      rel[boundary="administrative"][admin_level="6"](area.searchArea) -> .regio;

      foreach.regio->.b
      {
        .b map_to_area -> .newarea;
        (
            nwr[place="neighbourhood"](area.newarea);
        )-> .neigh;
        foreach.neigh->.n
        {
            make out _row = n.set(lat()) + " " + n.set(lon()) + ", " +  b.set(t["name"]);
            out center;
        }
      }
    );


texts: I am trying to label a line which contains 431 segments.  I need to label each segment with a simple GIS_ID number.  When I label the segments, many go unlabeled.  Is there any way I can label all 431 segments?


texts:
I'm trying to set some labels and styles for some vector layers that are loaded programatically from a Spatialite database. I have the following code, which returns no error, and includes the changes in the .qgis file, but yet doesn't show on the canvas:

uri = QgsDataSourceURI()
uri.setDatabase(database_path)
uri.setDataSource('', table_name, 'geometry')
layer = QgsVectorLayer(uri.uri(), display_name, 'spatialite')
QgsMapLayerRegistry.instance().addMapLayer(layer)

layer.setCustomProperty("labeling/enabled", "true")
layer.setCustomProperty("labeling/fontFamily", "Arial")
layer.setCustomProperty("labeling/fontSize", "10")
layer.setCustomProperty("labeling/fieldName", "source_id")
layer.setCustomProperty("labeling/placement", "2")

Can someone please point me in the direction of a working example? Examples of setting colors/widths etc. would be even more appreciated.

Thanks in advance for any help. Dan


texts: Algorithm for clipping large volume of points to irregular geometries?

texts:
I am currently using a home-brewed-hacked-together method for clipping a large number of points (~500k, could be more or less depending on the day) by a circle with a fixed radius. The output is a list containing only the points within the circle

This is my current method using Python & OGR bindings:

plot_points = []
circle = #< Some sort of circle geom>
for point in point_array:
     thing = ogr.Geometry(ogr.wkbPoint)
     thing.AddPoint(point[1], point[2])
     if thing.Within(circle):
         plot_points.append(point[0])

This method works, but is very slow.

I would like to avoid adding another dependency, as my list of dependencies is growing ever-longer.

Is there some sort of sorting algorithm I could implement that would speed up this process that uses ogr?

I have been looking into implementing quad-trees somehow, but have been unable to find any good examples of implementation in this context.

The points are initially stored in a 2D numpy array. It would be useful from an efficiency standpoint to not have to convert them into geometries.

I would prefer to clip to irregular polygons, not just circles.

Points are currently in UTM.


texts: Adding popup to Leaflet but only on click?

texts:
As there a way to add a popup to a marker but only show it when it is clicked?

As it seems to be default to show it when one is loaded.  This is no good when markers are getting updated.

Currently using this

 mymarker = new L.Marker(e.latlng, {
                icon: yellowIcon,
               draggable: true
             });

             map.addLayer(mymarker);
             mymarker.bindPopup('
You are here ' + username + '
').openPopup();


texts: Return rows where distance is within X meters

texts:
So how do I write a query to get all rows within X meters of a certain lat/long coordinate?

I tried looking into the docs, but it's super confusing, ideally I just want to write

SELECT *
FROM   locdat
WHERE  ST_DWithin(ST_MakePoint(-71.1043443253471, 42.3150676015829),2000, true )

but I don't understand anything of the syntax they use in the docs or something like select ST_MakePoint(loc_point) from locdat won't work either.

loc_point is data like (54.66965,20.64888) of type text.


texts: Exporting GeoJSON file with correct order of coordinates?

texts:
I was testing rendering of data via GitHub in geojson format, because I wanted to use it for INSPIRE data. INSPIRE data are in GML 3.2.1 format. I've downloaded one of datasets from http://services.cuzk.cz/gml/inspire/cp/epsg-4258/ (which is in ETRS). I needed to get json file from it, so I've opened GML file in QGIS (version 1.9) and saved it as GeoJSON file (CRS=EPSG::4326) and then uploaded to my GitHub. Order of coordinates in geojson is (easting, northing), but after saving file from QGIS it's (northing, easting). My data comes from Czech Republic, but it's rendered in Yemen.

Does anybody know how to switch order of coordinates (or axis) in GeoJSON file?


texts: Landsat8 NDVI values not in the -1 to 1 range

texts:
I was trying to create an NDVI from Landsat8 OLI image from July 2013.  I ran DOS1 and TOA.  The values range up to 3.  The original files were downloaded from EarthExplorer(EE), ran on QGIS and Multispec.  The original data was downloaded again from EE and then ran only in QGIS (band calculator).

What am I doing wrong?

I switched to a landsat 5 image and i am having some problem. After the TOA and DOS1, these are the results for visible Red:'RT_LT50010862008320COA00_B3.sta' Number of Bad Lines = 0 Total Number of Pixels = 2,969,703 Histogram Summary Table First Last Line First Last Column Line Line Interval Column Column Interval 1 7251 20 1 8181 1 Channel Data Range Mean Median Standard Description Deviation 1 0.0075 to 0.6560 0.0409 0.0432 0.0294


texts: Overlaying misaligned DGN file in ArcGIS for Desktop?

texts:
I received a DGN file, I believe created in Bentley 8, and there is a bit of a misalignment problem. The dgn, when overlain on a Bing Maps aerial or NAIP image, shows up about 20 m or so off. I'm set up in NAD83 State Plane 1601 Kentucky North and using ArcGIS 10.2.1 for Desktop.

It's important that these line up fairly close because this dgn will be used to create a regulatory boundary for field surveys.

So my question is, are the aerials correct and the DGN is just slightly off or vice versa? I assume the DGN is just slightly off. My initial inclination was to align the buildings in the dgn with the buildings in the aerial photograph and go from there. However, I started getting a little nervous that maybe that would end up putting one of the field workers off the mark.

Any suggestions?


texts: Creating polygons based on point layer using QGIS

texts:
I have a point layer of residential coordinates which includes an attribute field indicating the most visited hospital for the people living at the coordinate. From these data, I want to create a polygon layer of "hospital areas" that groups neighboring coordinates with the same most visited hospital, with a bit of an outside buffer, as illustrated below. Ideally, I want the polygons to cover also "weird" coordinate pairs, i.e. any coordinates which are clearly geographically located within a larger cluster with the same most visited hospital, but which they themselves have a deviating most visited hospital. Please see the example of Hospital B in the middle of the Hospital A cluster in my illustration.

How do I do this? The real data covers many thousand coordinate pairs and a total of 150 different hospitals.


texts: Counting points within radius of another point when a condition is met

texts:
I'm trying to create what I think would be called a spatial join in Qgis as follows:

I have two layers of points both with date fields

What I am trying to achieve is to create a count for every point in layer 1 of the number of points in layer 2 within a set radius (say 500 metres) of point 1.

Then I would like this count to be stored as a field in the attribute table for layer 1.

The issue I have is that I need this to occur when the two date fields match up i.e. when the observation in layer 1 has the date field 2005m1 then I only want to count the points in layer 2 when the date field in layer 2 matches up

My dates are stretched over a period of about 10 years so if there was some way to create a loop that automatically went through each subsequent month that would be brilliant.

I think this post references a similar process but doesn't allow for the join to be based on the additional date requirement:

Count points within each buffer QGIS

They mention that it could be achieved without buffering the points in layer 1 and this would probably be preferred but if anybody knows how to do it with buffers I'm equally happy.


texts:
I have two carstosat stereo pair images. All i need to do is to extract DEM from those.
So, i tried it doing with OSGeo4W. but i didn't succeeded in it.

So IS there any other softwares or opensource for extracting DEM from cartosat pair?

Thanks in advance


texts: Export to DXF from PostgreSQL using ogr2ogr

texts:
I have a large dataset in a PostGIS database (OS MasterMap). My aim is to export specific areas of this dataset to DXF, while retaining the attributes of the data as layers in the DXF.

My current ogr2ogr command looks like this:

$ ogr2ogr -f "DXF" os_mm_export.dxf PG:"host=11.11.11.11 dbname=db user=user password=xxxxx" -sql "SELECT ST_INTERSECTION(a.the_geom, (SELECT ST_SETSRID(geom, 27700) FROM schema.site_boundary)) AS geom FROM tbl_os_mmap AS a WHERE ST_INTERSECTS(ST_SETSRID(a.the_geom, 27700), (SELECT ST_SETSRID(geom, 27700) FROM schema.site_boundary))"

This works fine to export the data to DXF from the specified area, but the DXF only comes out with a single layer '0'.

What I want is to specify a 'line type' field to be used to create layers in the DXF.

There are also points to be exported, for those tables I would like to be able to specify which field would be used to label the points.

I realise it would be possible to break the process down into multiple steps, but I'm interested to see if it would be possible to do both things in a single ogr2ogr command, as this would obviously be the fastest and most elegant solution.


texts: Uploading 50Gb GIS data (csv files) contained in nested directories into PostGIS database

texts:
I have two years worth of GIS data contained in CSV files stored in nested directories (per year, per month) so it looks like:
2013
   01
       2013-01-01.csv
       2013-01-02.csv
       etc.
   02
   etc.
2014
   01
       2014-01-01.csv
       2014-01-02.csv
       etc.
   02
   etc.

The files contain the fields: ID, timestamp, lat, long. Combined, it is around 50GB worth of data. I am looking for the most efficient way to load data stored in a folder hierarchy into a single table in PostGIS. (This is my 'training data', and I'll be working with 500GB later on.)
At the moment I am thinking about writing a PyQGIS script based on os.walk (still have to work out how, tips are very welcome) for loading all the points into QGIS, saving them into the same layer, and then loading that layer into PostGIS. But I can imagine this would not be most efficient.
Any suggestions for alternative routes?


texts: Using SLD in GeoServer with OpenLayers 3?

texts:
I am facing display problem in SLD. Working on GeoServer 2.8 and OpenLayers 3.

I have one polygon layer in which I am applying SLD, when I view it in my web application polygon number is repeating each times and when I use centroid function it resolves.

But when I update SLD then that "centroid code" part vanishes from entire code.

attached screenshot with this (repeated questions highlighted in red circle)

searching on Google also

code link is: https://jsfiddle.net/5f84nwwj/

Update:updating SLD using Geoexplorer

please check the screenshot


texts: Show labels of a KML file in Google Earth

texts: I have a KML file created from an Excel spreadsheet which contains name, lat and long. How can I display the name in Google Earth?


texts: Reprojection causes polygon to stretch across globe

texts:
I've run into this issue a few times before in the past and I'd finally like to get a better idea of what's actually going on.

When reprojecting a layer, some polygons that are located around edge of the map extent will get stretched completely across the map. Is this indeed because of their location on the edge of the extent, and does this dramatically affect the accuracy of the polygon? Or is it more of a rendering artifact that does not seriously affect the underlying geometry?

The pink coral data in EPSG 4326 is the original data:

The reprojected (to equal area) orange coral data in EPSG 3410 is here:

Edit: Data w Pacific central meridian
It does appear to be a 180th meridian issue.


texts: Connecting to File Geodatabases from ArcGIS Pro?

texts:
I have File Geodatabases that I created in ArcCatalog and ArcMap 10.5.  I am trying to open maps in ArcGIS Pro from these File Geodatabases, but when using the "import Map" tool only folders show up, and no Geodatabases.  If I use the "+Connections >> Add Database >> Select Existing Geodatabase" tool my file Geodatabases show up and I can connect to them, but none of the .mxd files are shown, only feature classes, toolboxes, python scripts, and .csv's.

How do I connect to a .mxd file within a file Geodatabase in ArcGIS Pro?


texts:
After launching Qgis 3.2 Desktop I receive the following error messages:

    Impossibile caricare PyQt. Il supporto Python verrà disabilitato.

  Traceback (most recent call last):   File "", line 1, in  ImportError:
  DLL load failed: Impossibile trovare la procedura specificata.

  Versione Python:
  3.6.0 (v3.6.0:41df79263a11, Dec 23 2016, 08:06:12) [MSC v.1900 64 bit (AMD64)]

  Versione di QGIS:
  3.2.0-Bonn 'Bonn', bc43194061

  Percorso Python: ['C:/OSGEO4~1/apps/qgis/./python',
  'C:/Users/Mascia/AppData/Roaming/QGIS/QGIS3\profiles\default/python',
  'C:/Users/Mascia/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins',
  'C:/OSGEO4~1/apps/qgis/./python/plugins',
  'C:\OSGeo4W64\bin\python36.zip',
  'C:\OSGEO4~1\apps\Python36\DLLs',
  'C:\OSGEO4~1\apps\Python36\lib', 'C:\OSGeo4W64\bin',
  'C:\Users\Mascia\AppData\Roaming\Python\Python36\site-packages',
  'C:\OSGEO4~1\apps\Python36',
  'C:\OSGEO4~1\apps\Python36\lib\site-packages']

    Blockquote

and this

  Si è verificato un errore durante l'esecuzione delle istruzioni
  seguenti: qgis.utils.uninstallErrorHook()

  Traceback (most recent call last):   File "", line 1, in  NameError:
  name 'qgis' is not defined

  Versione Python:
  3.6.0 (v3.6.0:41df79263a11, Dec 23 2016, 08:06:12) [MSC v.1900 64 bit (AMD64)]

  Versione di QGIS:
  3.2.0-Bonn 'Bonn', bc43194061

  Percorso Python: ['C:/OSGEO4~1/apps/qgis/./python',
  'C:/Users/Mascia/AppData/Roaming/QGIS/QGIS3\profiles\default/python',
  'C:/Users/Mascia/AppData/Roaming/QGIS/QGIS3\profiles\default/python/plugins',
  'C:/OSGEO4~1/apps/qgis/./python/plugins',
  'C:\OSGeo4W64\bin\python36.zip',
  'C:\OSGEO4~1\apps\Python36\DLLs',
  'C:\OSGEO4~1\apps\Python36\lib', 'C:\OSGeo4W64\bin',
  'C:\Users\Mascia\AppData\Roaming\Python\Python36\site-packages',
  'C:\OSGEO4~1\apps\Python36',
  'C:\OSGEO4~1\apps\Python36\lib\site-packages']


texts:
I have around 1000 tiff files (each file around 300 MB in size), I published them via Geoserver as a mosaic using Image mosaic plugin; when trying to preview the published layer using openlayers, the layer performance is extremely slow at small scales (1:500000 and 1:1000000), and the server's CPU and memory are significantly consumed at these scales.
Accordingly, I am trying to cache this service using GeoWebcache that is shipped with Geoserver, I started the caching process using certain defined gridsets (2000000,1000000, 500000, 250000, 125000, 64000, 32000, 16000, 8000, 4000, 2000, 1000, 500). The problem is that the caching process has consumed all the space on my hard disk (e.g. the cache folder for 1:000 scale is larger than 400 GB) considering the fact that the whole mosaic size is less than 300 GB (with overviews).

The caching process is taking much much more space than the original data itself, so is this a normal behavior with GeoWebcache (I tried caching the same data at the same scales with ArcGIS server and it consumed less than 250 GB)?


texts: How to find the original vertex location from the OnVertexMoved event

texts:
So I am creating a fairly small, and hopefully simple, editor extension. This is in Arc 10.4.

What I would like to do is this. If a user is editing a feature and is manipulating the edit sketch, when they move a point I want any other coincident vertices in other features to move with it. Executing this should be easy enough:

Loop through all feature classes in the GDB
Buffer the original point and execute a spatial query
See if the original point exists on any discovered features
If a feature has the original point, move to the new location

#2 is the problem. The issue is that the OnVertexMoved event only provides the new point location, not the old. My first try was to get the selected feature from the Editor and then pull other information from the IEditSketch interface like so:

            //from inside OnVertexMoved
            IFeature source_feature = ArcMap.Editor.EditSelection.Next();
            IEditSketch e_sketch = ArcMap.Editor as IEditSketch;
            IGeometryCollection gc = source_feature.Shape as IGeometryCollection;
            IGeometry geo_part = gc.Geometry[e_sketch.Part];
            IPoint orig_point = (geo_part as IPointCollection).Point[e_sketch.Vertex];

But e_sketch.Part and e_sketch.Vertex are both invalid (-1). Is there a good way to grab the original vertex location? At this point the original feature's geometry should be unmodified so even just the part and/or index of the vertex should be sufficient.


texts: What is most accurate Spatial reference system for calculate Area entire world

texts:
I want to know what Spatial reference system is most accurate to calculate area?
I want to use it for entire world (-180,+180,-90,+90 area in WGS84 system) I have every thing ok I only need to reproject my layer to an accurate SRS and calculate area,Can you help me to find most accurate Spatial reference system for calculate Area entire world?


texts: approximation position in RINEX file header

texts:
In RINEX observation file header, there is a line with APPROX POSITION XYZ label.
Value of this line is approximate marker position in geocentric coordinate-system (unit: meter). Following Example shows a RINEX OBS File header. How
much is the standard deviation of These values?

     3.02           OBSERVATION DATA    M (MIXED)           RINEX VERSION / TYPE
TPP 3.10.4                              20180325 235959 UTC PGM / RUN BY / DATE
CKSV                                                        MARKER NAME
 -2956619.1390  5075902.2771  2476625.5647                  APPROX POSITION XYZ


texts: How to programmatically determine the size of a feature class in a file geodatabase?

texts:
My answer to the question How to get size of file geodatabase feature class on disk? shows how one can enable a column in ArcCatalog showing the file size of FGDB feature classes in the contents view.

Does anyone know how ArcCatalog is able to calculate this, and is there a method or property exposed in ArcObjects or another API to determine this programmatically?

I had a look through the esriGeodatabase ArcObjects namespace and did not find anything relevant, nor in the file geodatabase API, but I may have missed something.


texts: Exported shapefile from ArcGIS won't open in QGIS

texts:
I exported a shapefile from ArcGIS and wanted to load it into QGIS. However, the following error occurs:

Invalid Data Source: /Users/Desktop/Documents/merged_buffer.shp is not a valid or recognized data source.

When I check the files in the folder, the icons to the left seem white. The other shapefiles have a black icon, originate from QGIS itself and can be loaded correctly. What is the problem here?


texts: qgis set layer "unselectable" - Avoid identifying features in multiple layers

texts:
In MapInfo, it is possible to turn 'selectability' off on layers, which means that when you click on the map objects on that layer aren't selected. It's very useful.

Is there a way to accomplish this in QGIS?

I should add that what I am specifically referring to the information tool. In the attached screenshot you can see that the street, the neighborhood and the city were all highlighted, but I would like to restrict it to selecting roads only.


texts: Converting adjacent singlepart polygons to multipart?

texts:
I am interested in converting a polygon featureclass created via Create Fishnet (Data Management) from singleparts to multipart.  QGIS has a very handy tool called singleparts to multipart that does the job, however, I need to stay within the arcpy/Python realm.  The ArcGIS equivalent appears to be Dissolve, yet the tool does not preserve individual feature geometry when the polygons share a border.

The following is an example of the type of fishnet polygon data I am working with, which was generated from the below script.

import arcpy

# Create a fishnet with 9 columns and 9 rows
# with origin at (1, 1) and output geometry is set to default (POLYGON)
arcpy.CreateFishnet_management("C:/temp/fishnet2.shp", "1 1", "1 9", "1", "1", "9",     "9", "#", "NO_LABELS", geometry_type = "POLYGON")

How can I convert a fishnet polygon from singleparts to multipart using a scripting approach?

*Note. There is a similar question at GIS SE (Create a single connected polygon from multiple polygons), although it appears to be limited to desktop tools.  I am looking for any Python or arcpy-based approaches.


texts: Why is zero length polyline resulting from arcpy.da.InsertCursor()?

texts:
Below is a code snippet that I have extracted from something that I am working on.  It should create two polyline features: the first is expected to be tiny, while the second is expected to be small.

The small polyline is created as expected, but the tiny line is created as a feature of zero length when I think it should be about 4 metres long.

I have experimented with setFalseOriginAndUnits on the SpatialReference object to see if I can get a greater than zero length but have been unsuccessful.  I thought that this was worth trying on the basis of my learnings at Controlling XY Resolution on file geodatabase feature class output from Append?

How can I get this tiny feature to be created with a tiny but not zero length?

import arcpy

if arcpy.Exists("C:/Temp/test.gdb"):
    arcpy.Delete_management("C:/Temp/test.gdb")

arcpy.CreateFileGDB_management("C:/Temp","test")
sr = arcpy.SpatialReference("Geocentric Datum of Australia 1994")
sr.setFalseOriginAndUnits(-400.0,-400.0,10000000.0)
arcpy.CreateFeatureDataset_management("C:/Temp/test.gdb","testFD",sr)
arcpy.CreateFeatureclass_management("C:/Temp/test.gdb/testFD","testLines","POLYLINE")
cursor = arcpy.da.InsertCursor("C:/Temp/test.gdb/testFD/testLines",["SHAPE@"])
# Tiny line but NOT zero length
point = arcpy.Point()
array = arcpy.Array()
point.X = 153.06436725
point.Y = -27.49210825
array.add(point)
point.X = 153.06437825
point.Y = -27.49208025
array.add(point)
polyline = arcpy.Polyline(array)
cursor.insertRow([polyline])
# Small line
point = arcpy.Point()
array = arcpy.Array()
point.X = 153.0640
point.Y = -27.4920
array.add(point)
point.X = 153.0650
point.Y = -27.4930
array.add(point)
polyline = arcpy.Polyline(array)
cursor.insertRow([polyline])
del cursor

with arcpy.da.SearchCursor("C:/Temp/test.gdb/testFD/testLines",["Shape_Length"]) as cursor:
    for row in cursor:
        print "Length of line is {0}".format(row[0])

The output that results from the above code is:

>>>
Length of line is 0.0
Length of line is 0.00141421356242
>>>


texts: How to resolve "Not all features are shown" in ArcGIS JavaScript API

texts:
I need use a feature layer which includes 7000 features in SceneView. But the HTML page only shows part of data, and the error shows:

[esri.layers.graphics.controllers.SnapshotController] Feature limit exceeded on layer "Community". Not all features are shown.

The MaxRecordCount is set as 10000.


texts: If Then Statement, Python Label expression, ArcGIS Pro

texts:
I'm trying to create a Python label expression in ArcGIS Pro for a parcel layer where I only want to label the parcel when a street number in the Address field is present.If a street number is not present I want label with a statement to identify missing address. I am labeling the address of the parcel. The field type is string.

Presently I was able to create an expression where I have identified NULLS and spaces in the field(see below) but I'd like to do the same for numbers since I have fields where street numbers are not present but streets names are present.

def FindLabel ( [ADDR],[ACRES_GIS] ):
 If ([ADDR]) in (None, " ","","  ") :
  return "<BOL><CLR red='200'>"+'ADDR Absent'+ "</CLR></BOL>" +"\n"+"ac: "
  [ACRES_GIS]
 else:
  return "Addr.: "+[ADDR].split(" ")[0]+"\n"+"ac.: "+[ACRES_GIS]

Any suggestions on how to format such a script as the one I am describing?


texts: Publish a Feature class on ArcDesktop with C#?

texts: I have a feature class in ArcCatalog.  I need to save the feature class to a Map and then publish the map to the Arcserver.  All programatically and hopefully in C#.  I can't seem to find any examples at ESRI or by googling.  Any help and examples of this would be absolutely amazing. But just leading me down the right path is cool too.


texts: Sources for meterological data with high resolution

texts:
I was running Maxent on an endemic species of andaman and nicobar.But the climate data
I have from worlclim.org ( tile number 29 ) is not of high resolution and thus I can't see the distribution clearly.Can anyone please suggest me where i could get environment data for
andaman and nicobar islands of high resolution.
The middle left hand side with purple dots is the andaman nicobar.

The current resolution is 1km and i want better one like 10 m or so , so that I can see small landmasses also with detail and clarity.


texts:
My company has workers that use Windows and Linux to manage a QGIS project. However, when we use the same project created in different operate systems the layers lose file paths.

Is there possibility to create a compatible project created in Linux or Windows without the necessity to reset layers' file path?

Example below:

But in Linux the path of first one layer Americasul_x... is /mnt/mairi/GEO_DADOS/3_VETOR/3.1_SHP/BASES_GERAIS_SHP/Americasul_xxxx_xx.shp and it is different of how Windows system set it.


texts: Adding feature to highlight layer forces redraw of original layer

texts:
I have two VectorLayers setup, one for my actual features and another one for highlights (I'm rendering ~10000 features so I have a separate layer for highlight as recommended by ahocevar).

However, whenever I add/remove features to the highlight layer, it seems that the original layer is forced to be redrawn again, even though there was no interaction with that particular layer (besides hit detection).

Am I doing something wrong or does OL always rerenders all the layers when one of them changes?


texts: Create own "toolbox" or "toolset" in QGIS and move tools to that toolset

texts:
In ArcGIS Tools are located in "toolsets" wich are located in "toolboxes". One can create a new Toolbox and toolset and copy Tools from other toolsets into to the new toolset.

Is there a way to do the same in QGIS? To create a new "Provider" and copy Tools into that Provider?


texts: Prevent user from drawing a feature with Leaflet.draw

texts:
In my application using Leaflet.draw the user can set point markers on a map.

I would like to allow the setting only in a certain area.

The function isMarkerInsidePolygon(marker, polygon) returns true, if the marker is inside the polygone, false if outside.

For this purpose, I think I have to check after draw:drawstart, but before draw:created, if isMarkerInsidePolygon(marker, polygon) returns true or false.

If it returns false, I want to prevent setting the marker, but still keep the drawing mode active. Is there a posibility to prevent the setting of a marker in leaflet.draw?

There is a similar question asked by Ben, but no exactly the same case:
 Leaflet draw: prevent drawing through a custom overlay


texts: OpenLayers - changing layer type (wms to OSM) and WKT transformation problems

texts:
I've created the very simple tool which utilizes OpenLayers.Format.WKT with DrawFeature - polygon and ModifyFeature together, so I can create and modify, load and save polygons in/from WKT format.
I've had it working on OpenLayers.Layer.WMS, but now I'd like to work on OSM map as I need the city names to be shown.

The problem is that after changing the layer type, the polygon is loaded properly - it's displayed somewhere else than it is supposed to be. Here is my loading code (simplified)

var wmsLayer = new OpenLayers.Layer.WMS( "OpenLayers WMS","http://vmap0.tiles.osgeo.org/wms/vmap0?", {layers: 'basic'});

// when using this, it doesnt work properly: var wmsLayer = new OpenLayers.Layer.OSM( "OpenLayers OSM");

var features = wkt.read(data);
features.geometry.transform(map.displayProjection, map.getProjectionObject())
vectors.removeAllFeatures();
vectors.addFeatures(features);

I believe this is due to transformation problem, but I don't know how to fix it...


texts: Missing argument in SGS simulation of HPGL

texts:
Here is a simple example to test the Sequential Gaussian Simulation in open geostatistics library HPGL. I use the following file as hard data for the kriging:

NEW_TEST_PROP.INC

here is my simple code:

from geo_bsd import *
size = (55, 52, 1)
grid = SugarboxGrid(55, 52, 1)
prop = load_cont_property(NEW_TEST_PROP.inc, -99, size)
cov1 = CovarianceModel(type=1, ranges=(10,10,10), sill=0.4)
sgs_result = sgs_simulation(prop, grid,radiuses = (20,20,20),max_neighbours = 12,cov_model = cov1,seed=3439275)

However, I always get the following argument error:

Traceback (most recent call last):

File "", line 5, in

sgs_result = sgs_simulation(prop, grid,radiuses = (20,20,20),max_neighbours = 12,cov_model = cov1,seed=3439275)

File "/usr/lib/python2.7/dist-packages/geo_bsd/geo.py", line 272, in new_f
return f(*args, **kargs)

TypeError: sgs_simulation() takes at least 7 arguments (6 given)

I know it is trivial but I just cannot see what I'm doing wrong.


texts: Publishing a large and complex map to the web

texts:
I have created a complex map of my area and I wish to publish it as an interactive map on the web. The dataset appears to be too large for QGIS2Web to work - the script completes eventually but in loading into the browser it displays the first view but then fails to load the rest of the data, change views, or display the legend, and so it appears to hang. Is there any other plugin or application that I can use?
I would like to publish a map that enables users to zoom in and out, select and deselect layers, and click on the various map items. Much like QGIS itself really.


texts: Calculating area weighted average of attribute in one shapefile using a separate grid shapefile

texts:
I have a detailed polygon shapefile with streets, structures/houses, and open areas as land cover data in the attributes table, and each land cover type has a corresponding manning 'n' value associated with it in another column in the table. The manning 'n' value is in decimal format.

I also have a 20'x20' grid shapefile from another program (FLO-2D) with specific grid numbers in the attribute table.  I would like to overlay the grid on the land cover file to obtain a weighted manning 'n' value for each 20'x20' grid based on the percentage of land cover type that lies within each 20'x20' grid.  There are approximately 400,000 grid cells in the shapefile.  Both shapefiles have the same projection.

I'm using ArcMap 10.4 desktop, but do have the spatial analyst extension. The final product I need is the 20'x20' grid shapefile with a weighted average value of manning 'n' for each grid number in the attribute table.

I've tried the zonal statistics and the zonal statistics as table, but the results did not make sense.  Is this possible with the 10.4 desktop version and only the spatial analyst?


texts: Turning simple Excel file into CartoDB map?

texts:
I have Field 6 in my Excel sheet as a column of Zip Codes.

How do I get them to show up on the map using CartoDB?

My map continues to be blank.


texts: OpenRailwayMap not functioning with OpenLayers 2

texts:
I am trying for several days now to integrate a basic OpenRailwayMap into a project using OpenLayers 2 API but it's simply not working. I followed all the steps mentioned here and ended up with something like this:

http://codepen.io/payamcf/pen/XKPrvm

What am I doing wrong?


texts:
I've come across an annoying irritation. I have a .SHP converted from MapInfo which consists of a large number of features, a lot of which are made up of a whole heap of combined/multipart polygons. I'm doing a basic edit in deleting nodes. I get down to 2/3 nodes and then it won't let me delete them no matter what I do. The undo/redo indicates nothing is happening. Have tried removing and bringing up again. No issues until the last few nodes - seems to work ok prior, edit, select a number of ways, click node tool, click, click again, goes blue and doesn't do anything. It shivers but that's about it. Any suggestions?

UPDATE: Ahaaa...."After selecting vertices for deletion, click the Delete key. Note that you cannot use the mActionNodeTool Node Tool to delete a complete feature; QGIS will ensure it retains the minimum number of vertices for the feature type you are working on. To delete a complete feature use the mActionDeleteSelected Delete Selected tool."


texts: Automatic removal of an extra vertices lying on a straight line or polygonal segment using QGIS

texts: I am looking for a QGIS tool that would allow me to remove extra vertices on a line or segment of a polygon. By extra vertices I mean vertices that lie on a straight line. The other vertices must remain as they are.


texts:
I'm trying to use GDAL to read-in a large .pbf file containing OSM data for London. The problem is that ogrinfo returns an error message, Too many features have accumulated ... as shown in the reproducible example below, or outputs a file that seems to have no features (or -1 according to the reproducible example below!):

gdalinfo --version
# $ GDAL 2.4.2, released 2019/06/28

# get data
wget http://download.geofabrik.de/europe/great-britain/england/greater-london-latest.osm.pbf

ogrinfo -oo OGR_INTERLEAVED_READING=YES *pbf lines > x2
# $ Warning 6: driver OSM does not support open option OGR_INTERLEAVED_READING
# $ ERROR 1: Too many features have accumulated in points layer. Use OGR_INTERLEAVED_READING=YES mode
ogrinfo -oo INTERLEAVED_READING=YES *pbf lines > x3
wc x3
# $  27  65 740 x3 # contains no features

A couple of questions:

Can others reproduce this behaviour and is it expected?
What command can be used to make GDAL output the full dataset?

See below for the full output of the command, which seems right according to the documentation at https://gdal.org/drv_osm.html :

ogrinfo -oo INTERLEAVED_READING=YES *pbf lines
INFO: Open of `greater-london-latest.osm.pbf'

 using driver `OSM' successful.

Layer name: lines
Geometry: Line String
Feature Count: -1
Extent: (-0.511482, 51.285540) - (0.335437, 51.693440)
Layer SRS WKT:
GEOGCS["WGS 84",
    DATUM["WGS_1984",
        SPHEROID["WGS 84",6378137,298.257223563,
            AUTHORITY["EPSG","7030"]],
        AUTHORITY["EPSG","6326"]],
    PRIMEM["Greenwich",0,
        AUTHORITY["EPSG","8901"]],
    UNIT["degree",0.0174532925199433,
        AUTHORITY["EPSG","9122"]],
    AUTHORITY["EPSG","4326"]]
osm_id: String (0.0)
name: String (0.0)
highway: String (0.0)
waterway: String (0.0)
aerialway: String (0.0)
barrier: String (0.0)
man_made: String (0.0)
z_order: Integer (0.0)
other_tags: String (0.0)


texts: Adding pins (graphics) to ArcGIS map in WPF?

texts:
I have a wpf app with an ArcGIS map, I used the extent property to zoom and pan to Canada (the location I want) and added the pins (images) by adding a graphics layer, see code below.

currently I'm using this map:
http://server.arcgisonline.com/ArcGIS/rest/services/ESRI_StreetMap_World_2D/MapServer

I'd like to use this instead:
http://server.arcgisonline.com/ArcGIS/rest/services/World_Topo_Map/MapServer

the problem is that the pins show up in the middle of the map, the 0,0 location, when I switch to the new map. Are graphics specific to a type of maps or am I adding graphics the wrong way? What's the best way to add images to map?

<Grid x:Name="LayoutRoot">
    <Grid.Resources>
        <esri:PictureMarkerSymbol x:Key="PushPin"
                                  OffsetX="11"
                                  OffsetY="39"
                                  Source="pushpin.png" />
    </Grid.Resources>

    <esri:Map x:Name="MyMap"
              Extent="-125, 40, -60, 60">
        <esri:Map.Layers>
            <esri:ArcGISTiledMapServiceLayer ID="StreetMapLayer"
                                             Url="http://server.arcgisonline.com/ArcGIS/rest/services/ESRI_StreetMap_World_2D/MapServer" />
            <esri:GraphicsLayer>
                <esri:GraphicsLayer.Graphics>
                    <!-- Halifax -->
                    <esri:Graphic Symbol="{StaticResource PushPin}">
                        <esri:MapPoint X="-63"
                                       Y="45" />
                    </esri:Graphic>

                    <!-- Toronto -->
                    <esri:Graphic Symbol="{StaticResource PushPin}">
                        <esri:MapPoint X="-79.5"
                                       Y="44" />
                    </esri:Graphic>

                    <!-- Vancouver -->
                    <esri:Graphic Symbol="{StaticResource PushPin}">
                        <esri:MapPoint X="-123"
                                       Y="49.5" />
                    </esri:Graphic>
                </esri:GraphicsLayer.Graphics>
            </esri:GraphicsLayer>
        </esri:Map.Layers>
    </esri:Map>
</Grid>


texts: How to enable “supportsAppend” in ArcGIS Online REST API

texts:
I am looking for some help using AGOL and the REST API.
We have a project where users are making new routes for a network. We need to be able to upload and append their GPX data to an AGOL feature layer (they will not be updating existing features). Because AGOL doesn't accept GPX we convert the GPX to a GeoJson file on the website server and then send this to the AGOL environment however we get an error 'Append not enabled' (sync and track changes are turned off), so how do we enable this?
I can do this exact process using AGOL so why does it fail using the REST API?
I have made sure that the field format is correct for the new GeoJson file by exporting the structure from the existing feature service and because we can do this process manually.
All help very welcome, we have been stuck on this for weeks and its key to the project.


texts: In QGIS, is it possible to save graduated styles to use on different columns?

texts:
I'm using QGIS to produce thematic maps of 9 quantities (given in 9 columns in the attribute table) which have values broadly within the same range and I would like to plot them using the same intervals using 'Graduated' style.

However, while I can save a style, it saves the column name with it. When I load the style to attempt to apply it to another column, as soon as I change the column, it resets the intervals, which is exactly what I would like it not to do. Is it possible to somehow keep the interval ranges while changing the column?


texts: Getting information out of an ArcGIS .lyr file

texts:
I'm using ArcGIS 10 and I downloaded a Raster layer from here: http://bioval.jrc.ec.europa.eu/products/glc2000/products.php

When I unzip it I get the following files:

1 Aux-File (15KB)
1 ArcGIS-Layer (14KB)
1 Tif Image (37MB)
1 .lyr file (XML Document) (6KB)
1 .tif file (XML Document) (15KB)
1 .tif.aux file (XML Document) (5KB)
1 .tif.ovr file (OVR File) (1.8MB)
1 .tif.vat file (DBF File) (2KB)

When I open the .lyr Layer in ArcGIS I get a nice colormap with a good legend (as can be seen in the Table of Contents). Looking at the Attribute Table I see only numbers, no description.

Where does the .lyr file gets its information from?! I need its information to create a legend and to know which number is what.


texts: How to implement a Gaussian function in field calculator?

texts: Are there any build-in functions like Gaussian or Epanechnikov functions in field calculator?  Or I have to write that by myself?


texts: Creating a mapset and PERMANENT location for multiple GeoTIFF files in GRASS GIS?

texts:
I have a 2 year MODIS time series, with 45 files, all in tiff format. However, I don't know how to create a mapset for all files. I can only do it for one file. How do I create for my entire time series in GRASS GIS?

I'm Using GRASS on Windows.

This is my tif files

I open grass gis and try to create a location for all files and click on the option read projection and datum terms a georeferenced datafile, but the option is one file at a time, I would like to create for all files at once.


texts:
If any of my ESRI REST map services are down I want to be able to tell the user, at minimum with an alert or message.

I thought I could use onLayerAddResult(layer, error) but it doesn't seem to be getting called if the layer based on a stopped/down service is added to the map, only when a layer is successfully added.

I see some examples using esri.request but it doesn't seem like I should have to do an extra call to the service since in the Firebug console I see these calls like this to each service:

  http:
  //myserver/ArcGIS/rest/services/Streets/Mask/MapServer?f=json&dpi=96&transparent=true&format=png8&layers=show%3A-1&callback=dojo.io.script.jsonp_dojoIoScript2._jsonpCallback

With a response of either:

  *dojo.io.script.jsonp_dojoIoScript2._jsonpCallback({"serviceDescription"
  :
  null,"mapName":"Layers","description":"","copyrightText":"","layers":[{"id":0,"name":"LMBCO_Mask","parentLayerId":-1,"defaultVisibility":true,"subLayerIds":null}],"tables":[],"spatialReference":{"wkid":2263},"singleFusedMapCache":false,"initialExtent":{"xmin":806353.163236593,"ymin":139639.306949227,"xmax":1165380.94101437,"ymax":327833.751393672,"spatialReference":{"wkid":2263}},"fullExtent":{"xmin":824755.941014371,"ymin":38597.6402825604,"xmax":1146978.16323659,"ymax":323319.862504783,"spatialReference":{"wkid":2263}},"units":"esriFeet","supportedImageFormatTypes":"PNG32,PNG24,PNG,JPG,DIB,TIFF,EMF,PS,PDF,GIF,SVG,SVGZ,BMP","documentInfo":{"Title":"","Author":"","Comments":"","Subject":"","Category":"","Keywords":"","AntialiasingMode":"Normal","TextAntialiasingMode":"Force"},"capabilities":"Map,Query,Data"});*

if the service is up and running, or:

  *dojo.io.script.jsonp_dojoIoScript2._jsonpCallback({"error":{"code":400,"message":"Unable
  to complete operation.","details":["Service 'NYCStreets/LMBCO_Mask' of
  type 'MapServer' does not exist or is inaccessible."]}});*

if the service is not available.

Is there a way to handle this reponse that is already coming back?

Thanks for your help!

Chad Wilcomb