texts: OpenLayers - Cannot draw a circle using createRegularPolygon

texts:
I have been trying this for a while but I cannot manage to draw a circle on a vector layer although my polygons are being created without any problems.

this.pj_epsg_900913 = new OpenLayers.Projection("EPSG:900913");
this.pj_epsg_4326 = new OpenLayers.Projection("EPSG:4326");

var bounds = new OpenLayers.Bounds(
        -13080498.908163622, 3833552.5668887002,
        -12978234.962928545, 3956763.6368140383
);

this.map = new OpenLayers.Map(canvas, {
    controls : [ new OpenLayers.Control.LayerSwitcher(),
                 new OpenLayers.Control.Navigation(),
                 new OpenLayers.Control.ZoomPanel()],
    maxExtent : bounds,
    maxResolution : 481.2932418958517,//156543.0399,
    numZoomLevels : 19,
    units : 'm',
    projection : this.pj_epsg_900913,
    displayProjection : this.pj_epsg_4326
});

this.osmLayer = new OpenLayers.Layer.OSM("OSM", null, {
    "tileOptions": {
        "crossOriginKeyword": null
    }
});

this.map.addLayer(this.osmLayer);

this.markersLayer = new OpenLayers.Layer.Vector("Markers", {
    styleMap : style_basic
});

this.geometryLayer = new OpenLayers.Layer.Vector("Geometry");
this.map.addLayer(this.markersLayer);
this.map.addLayer(this.geometryLayer);

Now so far so good. I then draw my polygon here and it works. I see my polygons I add to the geometryLayer. But when I try to add a Circle I can't see any.

{ "radius":"8.27781404016073", "centroid":"(-118.57075649707252,34.26082074825019)" }

var circle = OpenLayers.Geometry.Polygon.createRegularPolygon(
                new OpenLayers.Geometry.Point(lonlat_[0].substring(1), lonlat_[1].substring(1)).transform(this.pj_epsg_4326,this.pj_epsg_900913),
                circleMeta.radius.substring(1),
                30
        );

        var circleFeature = new OpenLayers.Feature.Vector(circle);
        circleFeatures.push(circleFeature);
var layer = this.map.getLayersByName('Geometry')[0];
// cleanup features first
layer.removeAllFeatures();
layer.addFeatures(polygonFeatures, circleFeatures);

So, the polygons are added but not the Circles. I double checked if the lat long is empty etc but its all there - only there is no shape drawn on the map.

What am I doing wrong? Anyone?


texts: Cul-de-sac (curve) turns generation in flow for parallel lines using ArcGIS for Desktop?

texts:
I'm using ArcMap 10.2 and I need to generate a parallel lines within the polygon and connect them one to other with the curve, with the proper radius and return radius (the same as cul-de-sac tool propose).

I found A good solution for making a proper amount of parallel lines within the polygon here Drawing parallel lines inside polygons (Well Paths)?.

But still I want to use cul-de-sac tool (Advanced editing tool) in the flow to generate such turns thru all the parallel lines I create within the polygon with this script, but I couldn't find any python description how to use this.

In ArcToolbox also there just Cul-de-sac Mask tool is available.

Maybe anyone can recommend how to solve it? Maybe there is a chance to use some curve generating instruments?


texts: value map setting" of a shapefile will be lost when you upload to a POSTGIS database?

texts:
I am busy making a skeleton shapefile, with many value maps either preset (in the shapefile itself) or either created by a link to a table.

I have also set up a POSTGIS Database where the shapefile and those extra tables will be stored.

Whoever I am concerend that when I am gonna upload this shapefile it will both lose the value "preset" value maps and the  value map relation link to my tables.

The reason I believe this is that when I had saved my skeleton shapefile as a different file to test something. I lost the value map settings.

Will this also happen when I upload the shapefile's to my POSTGIS database? If so I can warn them about they need to reset it manually everytime in Offline projects.


texts: Unique Natural Area Code (NAC) for a building

texts:
NAC is a geocode system for identifying an object anywhere on the Earth.

I am going to assign a unique NAC address  for each building in Detroit. But different NAC addressed can be assigned to each building, because the building is an area with a different range of latitude and longitude and so it can be assigned a range of NAC, not just a unique one.

Does any one knows how can I choose a single unique NAC for a building?

Or what software I can use to assign NAC addresses to buildings?

I have read somewhere that you can assign the NAC of the centroid of a building as a NAC of that building. But I do not know how much this is true.

(You can see a house here with different NAC addresses)


texts: Mapinfo to QGIS. How to separate features in the same TAB

texts:
Why does QGIS open my Mapinfo TAB and automatically merge all features (that causes  the colors of each features to change to a random unique color for all)?

How can I separate those features in QGIS and keep them in the same TAB like they are in Mapinfo... so I can apply the right color to each of them?


texts:
This is using ArcGIS10 arcpy with spatial analyst.

I need to call Viewshed, but with an observer height set. I have tried using a Z-aware point feature set, but that only throws an error. I am feeding an arcpy.PointGeometry object with a Raster object of a dem for my inputs.
e.g.

viewshed = arcpy.sa.Viewshed(dem, source_point)

I am using dynamically generated points as the source for my observer points for the viewshed. Feature classes have to be built in memory.

I have tried several ways to build an in-memory feature class with an OFFSETA field, but so far it has not worked.

Anyone have any idea how to use the Viewshed OFFSETA, or another way to set observer height when calling Viewshed inside python using an in memory feature class?

Here is the solution I came up with, using Matthew Snape's advice. dem, source_point, offseta, and offsetb are all inputs to function this is contained in:

   # Construct viewshed input
   arcpy.CreateFeatureclass_management("in_memory", "viewshedinput", "POINT")
   arcpy.AddField_management("in_memory\\viewshedinput", "OFFSETA", "FLOAT")
   arcpy.AddField_management("in_memory\\viewshedinput", "OFFSETB", "FLOAT")
   rows = arcpy.InsertCursor("in_memory\\viewshedinput")
   row = rows.newRow()
   row.shape = source_point
   row.OFFSETA = offseta
   row.OFFSETB = offsetb
   rows.insertRow(row)
   del rows, row
   # Calculate viewshed for source point
   viewshed = arcpy.sa.Viewshed(dem, "in_memory\\viewshedinput")


texts: OSX 10.9 + QGIS 2.0.1 - Missing osgeo module

texts:
Running the following:

QGIS 2.0.1 (via Kyngchaos)
OSX 10.9
GDAL/OGR 1.10.1 (via Homebrew).

Recently, I've become unable to use the OpenLayers plugin, which is complaining that I'm missing the osgeo module.

Has anyone else dealt with this?


texts: Error when using Azimuth and Distance plugin

texts: I've use Azimuth and Distance Plugin in the past and it worked with an earlier version of Qgis.  It even kind of worked this morning.  I realized I had the incorrect active layer.  I've closed the project and reinstalled the Plugin.  However, I still get this error.


texts: Scale-dependent displaying of a layer using OpenLayers plugin?

texts:
When I am creating a layer in QGIS the layer position/lines does not stay the same when I zoom in closer than, it seems like 1: 1000. I am adding 2 images to explain the problem. The moment I zoom in closer than 1: 1 000 then the original layer does not display correctly anymore.

Please, can you assist in this matter?


texts:
I have some grib files from ECMWF ERA-Interim Reanalysis in which all the data are projected onto the N128 Gaussian grid (256 points of latitude x 512 points of longitude, https://www.ecmwf.int/en/forecasts/documentation-and-support/gaussian_n128).

I need to compare these data with some that have been geo-located with a GPS and are therefore in WGS84 latitude and longitude.

How should I do?

Is there a reprojection to do?

If so, which one?

I am googling like crazy but I can't find any clear explanation.


texts:
GIS is brand new to me, but I've spent a couple of days researching so that I can at least ask somewhat intelligent questions.

I've decided to use QGIS
I need to do an interactive floor plan on the web and shapefiles fits the bill

Problem: I need to put the floor plan in as a backdrop (raster) to all the polygons that will represent the rooms. But, I don't know how to properly set the scale so that when I use the 'measure' tool it will actually come up with correct measurements.

Any advice?

Edit: Adding sample image to show what I'm doing. As you can see there's a floorplan that I've imported into QGIS as a raster layer. I'm going to be tracing the shapes of the rooms into a separate shapefile, which I then want to use on the web (using Silverlight) to interactively update the polygon fill colors representing temperatures in the various rooms. The problem I'm experiencing is that when I load the image separately into the Silverlight application and then load the shapefile a) The position and shape of the polygon looks nothing like it does in the GIS app b) It's in the wrong place altogether, but I think I might be able to solve that one myself. The shape that's coming out totally different is the worry?


texts: Creating a custom datum transform to convert from sphere to spheroid using ArcMap?

texts:
I have Lat/Lon data that is for a sphere, but I am trying to convert these points to an oblate spheroid. Currently the sphere-based data points are offset when projected onto a map that is based on the spheroid geographic coordinate system. I believe this requires a custom datum transform, but it is not included in ArcMap and I have to specify the transform.

Any idea on how to convert latitudes used on a sphere-based geographic coordinate system to one for an oblate spheroid?

The spheroids are for Mars. One is a sphere with radius 3396190, and one is a spheroid with equatorial radius 3396190 and polar radius 3376200. It requires a custom datum transform, but I am unsure how to define such a custom transform. How is this done?


texts: Reading d NetCDF4 format climate data observation in R?

texts:
I downloaded yearly Temp/Precip observation by the globe from here: Free open source climate data archive, and I download those for my research case: Terrestrial Air Temperature: 1900-2014 Gridded Monthly Time Series (V 4.01) and Terrestrial Precipitation: 1900-2014 Gridded Monthly Time Series (V 4.01). I need to clip out all data observation record of all German weather stations. However, I intend to extract average annual temperature and annual precipitation for each of these station coordinates of Germany.

But, when I tried to load those data and read them in R, I ran into problems. Respective data archive's website can be found look here. I think data format could be NetCDF4 format based on the description from above web link 4. I installed raster and ncdf4 package to read the data in R, but such operations were failed.

Here is how data look like:

list.files("stella/data/air_temp_1980_2014/", recursive = TRUE)

 [1] "air_temp.1980" "air_temp.1981" "air_temp.1982" "air_temp.1983"
 [5] "air_temp.1984" "air_temp.1985" "air_temp.1986" "air_temp.1987"
 [9] "air_temp.1988" "air_temp.1989" "air_temp.1990" "air_temp.1991"
[13] "air_temp.1992" "air_temp.1993" "air_temp.1994" "air_temp.1995"
[17] "air_temp.1996" "air_temp.1997" "air_temp.1998" "air_temp.1999"
[21] "air_temp.2000" "air_temp.2001" "air_temp.2002" "air_temp.2003"
[25] "air_temp.2004" "air_temp.2005" "air_temp.2006" "air_temp.2007"
[29] "air_temp.2008" "air_temp.2009" "air_temp.2010" "air_temp.2011"
[33] "air_temp.2012" "air_temp.2013" "air_temp.2014"

Here is what I did in R:

library(ncdf4)
library(raster)
list.files("stella/data/air_temp_1980_2014/", recursive = TRUE)

nc_open("stella/data/air_temp_1980_2014/air_temp.1980")
raster("stella/data/air_temp_1980_2014/air_temp.1980")

Here is error that raised by Rstudio:

> nc_open("stella/data/air_temp_1980_2014/air_temp.1980")
Error in R_nc4_open: NetCDF: Unknown file format
Error in nc_open("stella/data/air_temp_1980_2014/air_temp.1980") :
  Error in nc_open trying to open file stella/data/air_temp_1980_2014/air_temp.1980

> raster("stella/data/air_temp_1980_2014/air_temp.1980")
Error in .local(.Object, ...) : Couldn't determine X spacing

Error in .rasterObjectFromFile(x, band = band, objecttype = "RasterLayer",  :
  Cannot create a RasterLayer object from this file.

FYI, here I put the data on the fly for the test on your machine, files can be found on the fly, Here and Here.

How can I read those data in R?

Plus, I only need to pull out all station coordinates of germany from yearly Temp observation for the globe?

How can I make this happen in R?

Any workaround to do this?

Any more thoughts?

How can I read those data correctly in R (as data.frame object in R) ?

How to deal with this data?

Any thoughts?


texts: The ratio Pixel distance to geographical distance?

texts:
I have had this problem for several days and I have no idea to solve it.

The problem is that I am trying to find out the length and width of ships from a Bing map. I got the location point (lon, lat),the image zoom level and also the image size. The following image is an example. 800x800 image at (53.5211,10.0657)with zoomlevel 19.

After the image processing methods, I got the ellipse-shaped mask.

I then calculated the pixel distance of major axis = 472.86  and minor axis = 92.58 of the ellipse. The idea is to calculate the ratio of geographical distance divided by pixel distance and to try and get each pixel points and its corresponding geographical location by the conversion from given lat/lon in WGS84 datum to XY in Spherical Mercator EPSG:900913.

I now have two geographical point locations and calculate the distance (d_g) by using haversine formula. Also I have the pixel distance (d_p). Then the ratio is r = d_g/d_p. However,the detected length and width are not close to the real value.

I found the problem is in the conversion from given lat/lon in WGS84 datum to XY in Spherical Mercator EPSG:900913.

This is the function I used to convert.

function [x,y] = latLonToMeters(lon, lat )

% Converts given lat/lon in WGS84 Datum to XY in Spherical Mercator EPSG:900913"

originShift = 2 * pi * 6378137 / 2.0;

x = lon * originShift / 180;

y = log(tan((90 + lat) * pi / 360 )) / (pi / 180);

y = y * originShift / 180;

end

Is there another way to get the geographical location by given only center point (lon,lat), zoom level and size of image?


texts: ModelBuilder, Iterate by Field and Viewshed

texts:
I want to compute the independent Viewshed for 360 different points. If I run these all together, I have one viewshed output which shows how many times any given cell is visible, but this does not tell me WHICH Observer Point these cells are visible from. The Observer Points tool has a limit of only examining 16 points. Therefore, I want to run a single viewshed analysis for every one of my points, and then do post-processing to determine which Point can see which cell. However, these points cover a large area (120km x 30km) , and the viewshed tool is extremely slow. I need to restrict the viewshed for each point to only consider a pre-defined limited portion of the total area.

The data is currently stored as such:

Feature class containing all points. Each point has a unique identifier (Image)
A folder for each Image (the name of the folder being the Image number), within which is a DEM, which is labelled as DEM%Image%.

The Viewshed Tool needs two inputs: the raster to be analyzed and a point coverage.
I would like to be able to iterate BOTH of these:
A) iterate by Field so that each point is fed into the Viewshed tool by itself
B) iterate the raster that is used for each image.

I have been trying to use the "Iterate by Field" iterator plugged into the Viewshed tool, but it only lets me set it as a "precondition". Not entirely sure if I am approaching this the right way.

Update: I can get the Iterate Feature iterator to populate the "Input Point or Polyline Observer Features" part of the Viewshed Tool, but still trying to work out how to match a unique DEM to each input feature.


texts: Defining areas for inclusion in polygon conversion using ArcGIS for Desktop?

texts: Using ArcGIS 10.2 for Desktop, I am trying to create a polygon shape(s) from a raster image. Is there a way to indicate areas I do not want included in the shape. Similar to a buffer I suppose, but that will not do what I want it to do. I have the image below, top is an inset. I want the "main" areas to be within a polygon, but if a point is too far out it will be left out, this can be arbitrary depending on the location and position of other points. Is there any way to control this? I have looked at using Minimum Bounding Geometry, Convex hulls, etc. The main problem is I need to do this for over 200 files and aside from hand drawing a polygon around each shape, is there any better way to do it?


texts: I want to do SpatialLines programatically

texts:
This is my dataframe:

df <- data.frame(ID = rep(1:100, each = 2), x =runif(200, 10, 90), y = runif(200, 10, 90))

For each ID I need to create SpatialLines, I haven't had good luck.

I think I have to split the dataframe and then transform to Lines class and after that the SpatialLines using the splited dataframe. All that using inside of a apply function. Am I right? How would you do it?


texts:
I'm looking for a PostGIS/Postgres query thats a bit tricky.

I want to find all polygons within a specified radius and all neighboring polygons to the polygon where I have my coordinates

SELECT *
FROM polygons_wgs84
WHERE
ST_Intersects(geom,
 ST_Transform(
     ST_Buffer(
       ST_Transform(
         ST_SetSRID(ST_MakePoint(9.116447, 60.548203), 4326),
       3857),
     1000),
 4326)
)

And the polygons that ST_Touches the polygon that contains this point.


texts: Create points (shapefile) with ogr2og from an attribute table

texts:
I have a csv file containing lon,lat coordinates. I want to create a vector layer (points with the lon/lat coordinates) using ogr2ogr.

I did the following:

ogr2ogr -f "ESRI Shapefile" "vector.shp" "coords.csv" -oo X_POSSIBLE_NAMES=lon* -oo Y_POSSIBLE_NAMES=lat* -a_srs EPSG:4326 -nlt POINT -skipfailures

The result is a shapefile with no geometry values. Do I have to define anything else in the above command?

A sample of the coords.csv is as following:

Lon_rounde  Lat_rounde  ISO3    Year    GNFR_Secto  SourceType  CH4
-17.05  32.675  PRT 2015    A   A   57.69167458

I also tried to use the whole name of lon, lat (with the same result) :

> ogr2ogr -f "ESRI Shapefile" "vector.shp" "coords.csv" -oo
> X_POSSIBLE_NAMES= Lon_rounde -oo Y_POSSIBLE_NAMES=Lat_rounde -a_srs
> EPSG:4326 -nlt POINT -skipfailures


texts:
While running topology against two identical polygon (parcel) datasets in ArcGIS and QGIS, QGIS detects errors that ArcGIS does not. In four separate instances in my township-wide parcel dataset QGIS flags "invalid geometry" (screenshot 1) where two island polygons touch at a point within the larger containing polygon (screenshots 2 and 3). Notice that QGIS does not have an issue with a number of island polygons as long as they do not touch each other (adjacent parcel to the NE in screenshots 1 & 2).

Why is this an invalid geometry error?


texts:
Looking to make geo-reference pdfs (or tiff) of the sort that can be used by the Avenza application in R. Can anyone provide some advice on tools that may be able to convert an R plot into a georeferenced pdf or tiff?.

I have made a stab at my own function to do this (tiff), but while the resulting map is roughly in the right place, it is scaled slightly incorrectly.

map2tiff=function(tiff_out
                 ,dir_temp="c:\\temp"
                 ,ext #required, extent of first object plotted
                 #,crs=NA
                 ,res=300
                 ,xy_inches=6
){

  require("raster")

  #query graphics
  graph_in=recordPlot()
  par_in=par()
  dims_in=dev.size(units = c("in", "cm", "px"))
  mai_in=par_in$mai
  pin_in=par_in$pin
  plt_in=par_in$plt
  usr_in = par_in$usr
  din_in = par_in$din
  fin_in = par_in$fin

  #manually calculate mai - default mai does not account for legend when plotting a raster imager
  mai_in = c(plt_in[3] * fin_in[2],plt_in[1] * fin_in[1],fin_in[2] - plt_in[4] * fin_in[2],fin_in[1] - plt_in[2] * fin_in[1])*1.5

  #compute new graph dimensions etc.
  dx=xmax(ext) - xmin(ext)
  dy=ymax(ext) - ymin(ext)
  x_rat=min(1,dx/dy)
  y_rat=min(1,dy/dx)
  x_in = xy_inches[1]*x_rat + sum(mai_in[c(2,4)])
  y_in = xy_inches[1]*y_rat + sum(mai_in[c(1,3)])
  xrat= dx / x_in
  yrat= dy / y_in
  xyrat = min(xrat,yrat)

  #create new extent to include margins
  lx=xmin(ext) - mai_in[2] * xrat
  ly=ymin(ext) - mai_in[1] * yrat
  ux=xmax(ext) + mai_in[4] * xrat
  uy=ymax(ext) + mai_in[3] * yrat

  #export plot to a file
  tif_temp=file.path(dir_temp,gsub("[.]tif","_temp.tif",basename(tiff_out)))
  tiff(tif_temp,height=y_in,width=x_in,units="in",res=res)
  replayPlot(graph_in)
  dev.off()

  #load and georeference map
  r_map = brick(tif_temp)
  extent(r_map) = c(lx,ux,ly,uy)

  #write final geo-referenced raster
  writeRaster(r_map,filename=tiff_out,overwrite=T)
  unlink(tif_temp)

}

library(raster)

set.seed(55)

xy=data.frame(x=rnorm(50)*5000+200,y=rnorm(50)*15000+200)
coordinates(xy)=~x+y
rtest=raster(ext=extent(xy),res=50)
rtest[] = rnorm(length(rtest))
writeRaster(rtest,"c:\\temp\\tiff001.tif",overwrite=T)

plot(rtest)
plot(xy,add=T,pch=16)

map2tiff(tiff_out="c:\\temp\\maptiff001.tif",ext=extent(rtest))

rdat=raster("c:\\temp\\tiff001.tif")
rmap=raster("c:\\temp\\maptiff001.tif")

the rasters in the two plots should exactly overlap, but don't:

plotRGB(rmap)
plot(rdat,add=T,alpha=.5)


texts: Which suitable buffer distance to choose for a kernel density analysis to deal with edge effects?

texts:
I would like to run a Kernel density analysis on point features in ArcGIS 10.2 for hotspot detection, and to do it is recommended to apply a buffer zone to the study area in order to reduce the edge effects during the kernel analysis.
I would like to know how to do to figure out the most suitable buffer distance, knowing that my study area is about 331 000 km², my output cell size will be 130 m, and  my search radius will be about 150 km (distance of the first statistically significant peak z-scores given by an incremental spatial autocorrelation analysis).


texts: Aligning two non-coincident equi-resolution raster grids in ArcGIS Desktop?

texts:
I have two rasters: a digital elevation model (DEM) and a Landsat band. The grid resolution of both is 30 m, but the grids are not coincident (see image: the dark raster on top is the DEM; the light raster underneath is the Landsat band).

I want to perform some transformation of the Landsat band so that its grid aligns perfectly with the DEM grid. I've tried resampling the Landsat band to 30 m in ArcGIS using the DEM as a Snap Raster in the Environment settings, but this does not have any effect. I've also tried clipping the Landsat band to the extent of a vector polygon (that the DEM is also clipped to), while setting the DEM as the Snap Raster, but to no avail.

I suspect this is happening because both rasters are already at the same resolution (just on different grids), but I really need their grids to be aligned to be able to do any further analysis.

Can I do this in ArcGIS Desktop?


texts:
I am returning a JSON to a GeoExt app and this is the code that is collecting the return,

var control = new OpenLayers.Control.Click({
                     trigger: function(evt) {
                      var lonLat = evt.xy;
                      if (popup) popup.close();
                      var  request = Ext.Ajax.request({
                      url:"http://128.196.142.94:9000/info?lat="+ panel.map.getLonLatFromViewPortPx(lonLat).lat +"&lon=" + panel.map.getLonLatFromViewPortPx(lonLat).lon,
                      disableCaching: false,
                      success: function(response) {
                              var obj = Ext.decode(response.responseText);
                              console.log(obj);
                              console.dir(obj);
                              addToPopup(lonLat, response.responseText);
                          },
                      failure: function(response) {
                              console.log('server-side failure with status code ' + response.status);
                              addToPopup(lonLat, "Server side Failure");
                          }
                             //callback: myhandler
                 });
}});

I wish to obtain the returned JSON and then assign it to the popUp but whenever I try to connect to the App it is returning the following error,

OPTIONS http://128.196.142.94:9000/info?lat=55.44873046875&lon=-118.17626953125 405 (Method Not Allowed)
XMLHttpRequest cannot load http://128.196.142.94:9000/info?lat=55.44873046875&lon=-118.17626953125. Origin null is not allowed by Access-Control-Allow-Origin.
test.html:222server-side failure with status code 0

If you simply pass the url http://128.196.142.94:9000/info?lat=55.44873046875&lon=-118.17626953125 you will obtain the JSON but for some reason there is failure that keeps occurring in the application.

Is there CORS that can be used in the OpenLayers or GeoExt. I have tried using ExtJS JSON-P but I am not sure what should I enter in the 'callback' of,

<script src="http://domainB.com/users?callback=someCallback"></script>


texts:
I am reading the Proj4js User Guide but cannot quite grasp it. I import the library

<script src="lib/proj4js-combined.js"></script>

Then,

function init() {

    (...) // reading from GML string - one polygon in EPSG:3031
    var features = format.read(strGML);

    Proj4js.defs['EPSG:3031'] = '+proj=stere +lat_0=-90 +lat_ts=-71 +lon_0=0 +k=1 +x_0=0 +y_0=0 +ellps=WGS84 +datum=WGS84 +units=m +no_defs';

    var source = new OpenLayers.Projection('EPSG:3031');
    var wgs84 = new OpenLayers.Projection('WGS84');

    OpenLayers.Projection.transform(features[0], source, wgs84);

    (...) // show polygon on map

}

Well, this doesn't work. Reading of GML string works fine as the polygon shows correctly on map if the coordinates are in WGS84 (and without the transformation, of course). But when I put in the string with other coordinates and try to transform it doesn't work. Is it possible to transform a polygon at once or do I have to do it point-by-point (and how would I do it)? Thanks

EDIT:
My GML string loks like this:

<gml:featureMember xmlns:gml="http://www.opengis.net/gml xsi:schemaLocation="http://www.opengis.net/gml http://schemas.opengis.net/gml/3.1.1/profiles/gmlsfProfile/1.0.0/gmlsf.xsd" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
     <feature:feature xmlns:feature="http://example.com/feature">
          <feature:geometry>
               <gml:Polygon>
                     <gml:exterior>
                           <gml:LinearRing>
                                 <gml:posList>591674.39 5022898.05 545682.5 4722908.1 571701.44 5322909.29 651691.25 5022904.6 591674.39 5022898.05
                                 </gml:posList>
                            </gml:LinearRing>
                     </gml:exterior>
               </gml:Polygon>
          </feature:geometry>
     </feature:feature>
</gml:featureMember>

Should I, instead of

OpenLayers.Projection.transform(features[0], source, wgs84);

write

OpenLayers.Projection.transform(features[0].geometry, source, wgs84);

or even something else?


texts: Setting up MouseWheelZoom on the map in OpenLayers3

texts:
Well, I`m working on my own pet-project and using for this OL3. I need much help with creating a zoom for my map. The documentation for OL says that they do not have a built-in zoom by default, and an example of using a polyphile is also given. Unfortunately, the code was not working, and I could not solve the error.

Is there any other way to arrange zoom for the map?


texts: QGIS over WatchGuard Mobile VPN with SSL

texts:
something quite strange occured today. I connected to my office through WatchGuard Mobile VPN with SSL enabled. I got loads of errors about SSL certificates, WMS errors (there were WMS layers in the QGIS project file) and eventually, QGIS 2.01 crashed.

I got some colleagues to connect via VPN from one of our other offices and the same happend. QGIS is installed on my work laptop and this happens when I try to open figures on our network servers via VPN. When I am in the office, with the same laptop and same QGIS installation, and I try to open the same file, these errors do not occur. Therefore the problem is somewhere in the VPN connection. Why would QGIS fall over the VPN connection?

Thanks
Antonio


texts: How to convert an OSM relation to a list of lat/lng points?

texts:
I would like a convert a relation (e.g. https://openstreetmap.org/relation/2851613) to a list of coordinates that can be displayed on a map, potentially using geojson. I ported over this converter https://wiki.openstreetmap.org/wiki/Relations/Relations_to_GPX (the converted code is at https://github.com/e-mission/e-mission-eval-public-data/blob/master/spec_creation/autofill_eval_spec.py#L71) and it works iff the relations are sane.

However, many relations are malformed. Although the relation spec https://wiki.openstreetmap.org/wiki/Relation says that the members have to be ordered, this relation for the route of bus #52 (https://www.openstreetmap.org/relation/8094017) has the following way order

Way West El Camino Real (516422275) as forward
Way Castro Street (417101895) as forward
Way North El Monte Avenue (52869637) as forward

which is clearly incorrect since since the 52 goes from Castro to El Monte via El Camino. So if I use the algorithm above, I end up with a map like so

At the same time, OSM displays the route without any zigzags, so it is clearly able to order the relation members correctly even if they are not ordered originally.

In response to @neogeomat, I looked at the nodes displayed by overpass and the nodes displayed by the reconstruction algorithm above for a correctly created relation, and it looks like overpass returns fewer nodes than the algorithm above.

Does anybody know how OSM displays relations, or have a pointer to an algorithm that deals gracefully with malformed relations?


texts: Force ArcCatalog to refresh TreeView and View after programmatically creating a table

texts:
I'm programmatically creating a table in a geodatabase and am trying to force ArcCatalog to refresh the display to show the new table.  I've tried the following which does not work:

Dim gxApp As IGxApplication = TryCast(m_application, IGxApplication)
... 'Create table here

'Then trying to force tree view to show new table...
gxApp.Refresh(0)

'And trying to force view to display...
gxApp.View.Refresh()

Missing something simple here.


texts: How to hide labels in ArcGISDynamicMapServiceLayer

texts:
Does anyone know how I would go about hiding labels in a  ArcGISDynamicMapServiceLayer?

This is part of my code so far:

var routeLayer = ArcGISDynamicMapServiceLayer(...);
map.addLayer(routeLayer);
var drawingOptions = new LayerDrawingOptions();
drawingOptions.showLabels = false;
var options = [];
options[0] = drawingOptions;
routeLayer.setLayerDrawingOptions(options);

After doing this however the layer isn't shown on the map.

Do I need to set the renderer property for drawingOptions as well?


texts: Creating single feature to represent group of individual features, based on location and # using ArcGIS Desktop?

texts:
I am looking for an algorithmic way of handling this situation, as I can do it manually with the ruler and some time, but it would be nice to find a programmatic to have Arc do the work for me.

I have a bunch of point features, and these point features can be grouped into a new single feature, but only if A) the individual features are <20m apart and B) there are >4 individual features within that 20m of each other (not just 20m from a centroid).

I can buffer the individual features 20m to see which neighbours they touch, but then I need to differentiate those overlaps into 20m cluster groups.

What I don't want is a situation where PointA, B, C, D, E are all 20m from each other and should be grouped, but PointF is <20m from point E but >50m from PointA and there are no other points around, so I need a method which can understand when an individual point should be grouped because it is <20m from another point, but then figure out if there are enough other points around to make a group (>4> and then decide which group it should belong to, if there are two grouped clusters it could fit into, based on which centroid it is closest to.

Does that make sense?

I'm using ArcGIS Desktop 10.0.


texts: Replacing Z-field in vector file with attribute table column for export as Z-enabled KML using QGIS?

texts:
Using QGIS, I have been unable to take a height field contained within the attribute table of a vector file, and make that height field the z-value of that point.

QGIS doesn't seem to have the ability to create/alter Z-values within the core program, I would rather not have to deal with SpatiaLite databases if I don't have to, which is the only other method I have heard mentioned.

Also tried v.to.3d as suggested here
Features to 3D from 2D in QGIS

and got either an empty output, an OGR "Data source is invalid ()" message, or both.  I am working with a shapefile if that is relevant.

Any ideas?


texts: Determining unknown Coordinate System using QGIS?

texts:
I'm trying to figure out the coordinate system of some data I downloaded from a no-longer-maintained USDA website.

The data's still accurate, it's just floating in space.

I've tried several state-plane systems, a number of UTM systems, etc.  My brute force method is running out of steam.

I'm using Quantum GIS, so what I've done is started a new project, imported the state of new mexico's counties to the system and let the metadata assign its coordinate system (EPSG:4269 - NAD83).

I then enabled on-the-fly reprojection, set the project's coordinate system to "NAD83 / UTM zone 13N," and imported my un-defined coordinate data.

I tried about 30 different coordinate systems so far, and none of them get the data any closer than almost-in-line east-to-west and out of line by about 1000 miles north-to-south.

The image below is a screenshot of qGIS with Erodability due to wind data, with an assigned system (EPSG:26915 - NAD83 / UTM zone 13N) based on a guess, which is obviously wrong.

I've marked the origin with a small dot, if it helps.

For what it's worth, I went over the metadata files, and I didn't see anything relating to coordinates.


texts:
I want to know how to make an anaglyph of a DEM using QGIS. I've been trying with Orfeo Toolbox but I couldn't.

Can I make anaglyphs with other toolbox or plugin?


texts: Find polygons that does not contain (any) polygons with PostGIS

texts: New to PostGIS, I have a polygon layer of water catchment areas which are formed by the specific monitoring points along the river upstream. For example, on one river, there could be multiple areas, all overlaping with each other. Or there could be a small area not overlaping with any other areas. The task is to find the most upstream catchment areas (red), the ones that does not contain any other areas. I have tried several geometry expressions, but could not get the correct result.


texts: Overwrite the input feature with the output feature with ArcPy

texts:
Let's say that I have 4 features classe:

Layer1
Layer2
Layer3
Layer4

I want to erase my LayerX with all the LayerYwhere Y<X.
So for example the erase feature for the Layer3 will be Layer2 and Layer1.

I'm using the function Erase_analysis.

Erase_analysis (in_features, erase_features, out_feature_class)

I don't want to create tousand of new layers, so during my iteration I try to use the same input name and output name.

My code using arcpy:

arcpy.env.overwriteOutput = True

for x in n:
    for y in n[:(x-1)]:
        arcpy.Erase_analysis('layer'+str(x), 'layer'+str(y), 'layer'+str(x))

But, and it was expected, the python consol give me this error

  Output Feature Class is same as Input Features, Failed to execute
  (Erase).

Is there a solution using arcpy ?

Is there a more pythonic way to solve this problem ?


texts: Are polygons created in Arcpy built in any order?

texts:
My script uses a search cursor to evaluate polyline coordinates and build polygons based on those coordinates.  The script works well, but I am curious as to whether anyone knows the order (if there is any) that polygons are built.  I know that the output polygon feature class reads the rows (and their coordinates) in order, but I do not know if they build them in order.  I could see the polygons being built according to which order they were read, or based on size, etc.  Any input would be appreciated.

Thanks in advance.


texts: Openlayers2 checking coordinate system of WKT or KML file

texts:
In my Openlayers2 based project I give the possibility for users to upload their own features in KML file or with a WKT string.

I only want to accept features with ESPG:4326 coordinates. Otherwise a error message should be displayed. "Please use only 4326 coordinates bla bla bla"

The issue is that I can't find any way to identify the input coordinates system to make the necessary condition.

Is there a way to identify the EPSG only with the coordinates?


texts: Making a heightmap DEM for use in putting together a relief

texts:
I'm new to GIS, and I'm just doing this as a hobby.
I've seen 3d relief maps like this one: https://i.redd.it/peq57uoa28c51.jpg and I wanted to try my hand at making my own. I used the topoView service from USGS, and was able to acquire a .tif file, but that's the furthest I got.
Does anyone know how to create a relief like the one in the picture above? It doesn't look like the original author has any tutorials or guides.


texts: Why do my points end up in Africa? (should be Croatia)

texts: I have the correct coordinates in a excel file, but when I import the file, the point always goes to the wrong place. (I need it in Croatia, and it puts the point in Africa) I've tried chaging CRS system several times, but nothing seems to be working... It always goes somewhere else like the pacific or america... Can anyone help?


texts: How to properly utilize DWITHIN Spatial Filter in OpenLayers

texts:
Given the following code:

var points = [];
points.push(new OpenLayers.Geometry.Point(-110, 38));
points.push(new OpenLayers.Geometry.Point(-108, 40));
var lineString = new OpenLayers.Geometry.LineString(points);
var filteredPoint = new OpenLayers.Geometry.Point(-108.25, 39.78);
var newFilter = new OpenLayers.Filter.Spatial({
        type: OpenLayers.Filter.Spatial.DWITHIN,
        distanceUnits: 'm',
        distance: 92600,
        value: lineString
}); //92600 meters, or 50 nautical miles
alert(newFilter.evaluate(filteredPoint));

When the alert is fired, it results in false.  I'm a beginner to OpenLayers, looking for a way in Javascript to determine whether or not a point is within a certain distance, either side of a line.  I cannot use a bounding box for this feature, because I want the user to be able to pick any two points and a distance from the line, so an MBR isn't available here.  filteredPoint is less than 20nm from the 2nd point in the linestring, so my assumption is that this alert would return true, but alas it is false.
Meanwhile, given the following code, in addition to that above:

var bounds = new OpenLayers.Bounds();
bounds.extend(new OpenLayers.LonLat(-110, 38));
bounds.extend(new OpenLayers.LonLat(-108, 40));
bounds.toBBOX();
var newFilter2 = new OpenLayers.Filter.Spatial({
        type: OpenLayers.Filter.Spatial.BBOX,
        value: bounds
});
alert(newFilter2.evaluate(filteredPoint));

That code gives a result of true inside of this bounding box.
I would think that the two bits of code are trying to find points in the block below represented by the asterisk.  Am I wrong in interpreting DWITHIN as such?

----------                     /
|      * |                 *--/
|  bbox  |                   /
|        |                  /
----------                 /

Any tips or guidance in the right direction would be greatly appreciated.  My overall goal is to determine how many points out of ~3700 points I have fall within a distance of the line.


texts: What technique is used to acquire client location?

texts:
So far I know how to grab the client location to a application or web application by using the users geolocation or GPS location. I remember the possibility of finding the users approximate location by using their IP address. How is this done if you know how it is done can you link me to a tutorial or spec in addition to a short description?

Edit (Added Picture):

Edit: Possibly found through ISP?


texts: From coordinates find country using Python

texts:
I have tweets with coordinates. I'm trying to bucket the tweets by country using the coordinate information. I know that I can achieve this using some hosted web services but all of these have API limits that I am sure to run into.

UPDATE:

I started writing this code learning from another resource about osgeo library.

import osgeo.ogr
shapefile = osgeo.ogr.Open("geodata/worldborders/TM_WORLD_BORDERS-0.3.shp")
layer = shapefile.GetLayer(0)

for i in range(layer.GetFeatureCount()):
    feature = layer.GetFeature(i)
    countryCode = feature.GetField("ISO3")
    countryName = feature.GetField("NAME")
    print countryCode, countryName

    geometry = feature.GetGeometryRef()

The geometry variable has a member function called Contains that can take another Geometry and checks whether the other is contained inside itself.

I have a coordinate value from Twitter data in the following format: [45.29680725, -75.92681408]. However an example on osgeo showcases creation of a point like this:

point = osgeo.ogr.Geometry(osgeo.ogr.wkbPoint)
point.SetPoint(0, 474595, 4429281)

How do I convert my coordinates into the above required format? If I do manage to convert, is this the right approach of mapping twitter coordinates to countries?


texts:
I want to write cubic meters in QGIS 3, I searched a lot but I didn't find any instructions.

Is there a way to write cubic meters as m³?


texts:
I have a cql filter as SQL. Im trying to parse this cql to OpenLayers.Filter using OpenLayers.Format.CQL.read() method. It all works fine since my input is ascii. When I type Turkish characters as VALUE an error occurs. I modifed regex belong to VALUE property in CQL.js as /(^('.+'|\d+(.\d*)?|.\d+))/
Now it accepts non-ascii characters but generated Filter is wrong. The regex i modified might affect somewhere.

What i try and what message i get actually like follow

var format_cql = new OpenLayers.Format.CQL;
format_cql.read("NAME = 'ŞAKİR'");

Error: ERROR: In parsing: ['ŞAKİR'], expected one of: VALUE: /^('\w+'|\d+(.\d*)?|.\d+)/

If i type SAKIR instead of ŞAKİR it returns a OpenLayers Comparison object well.
It also doesnt allow to type space character in value.

Could someone please help me ?

Regards


texts: QGIS 3 - Polygon symbology RGB based on different columns

texts:
I want to symbolize several values from one layer using R, G, and B values. So, column 1 would be the blue aspect, from 0-255 based on min and max values, column 2 would be green, and column 3 would be red.

I can't find a way to do this in QGIS 3 without creating 16 million rules. If I make multiple layers with color and transparency, then whichever layer is added last dominates the color spectrum.


texts: Buffering multiple feature classes from dataset based on length using ArcPy?

texts:
I am a beginner at Python and for my current project I am expected to list all polyline feature classes within a dataset, find the length, and then perform a buffer analysis based on it's length. This is what I have currently:

import arcpy
from arcpy import env
env.workspace = "X:/Oil_Project.gdb"
#setting workspace
fc=arcpy.ListFeatureClasses("","Polyline","Well_Data")
#create variable to find polyline features within well data
for features in fc:
    fcdescribe=arcpy.Describe(features)
    print "Name: " + fcdescribe.name
    #prints name for each class
    for length in arcpy.da.SearchCursor(features, ["SHAPE@LENGTH"]):
        print("Length in feet = {0}".format(length[0]))
        #finds length
        if length < 3000:
            buffer_width= 75
        elif length < 4000:
            buffer_width= 300
        else:
            buffer_width= 175

I really have no idea how to perform the buffer analysis on each individual polyline feature automatically. I'm also not sure if I can do anything with the buffer_width variable.


texts: Copying shapefiles within windows explorer (or - extension)

texts:
Is there a user-friendly way to copy shapefiles to a different file location?

ArcCatalog is not fast enough for this purpose, creating  new map and navigating to this new workspace takes a lot of time..

Does an extension for windows explorer (version 7) exist  (maybe specifying an advanced filter) ?


texts:
Yesterday I asked what would be the best method of creating an overlay in a web map that would allow for some interactivity.
Create interactive map with no server
Thanks to a variety of people who recommended I go down the D3 and Topojson route.
The example that looks like the perfect solution is http://bl.ocks.org/mbertrand/5218300
So I created a topojson from an ESRI Shapefile using the following command
topojson --id-property full_name -o southeast.json southeast.shp

I have since tried every possible permutation and example to try and get the topojson loaded and I have hit a brick wall. I just do not understand what I am doing wrong.
I have used the following examples
http://www.schneidy.com/Tutorials/MapsTutorial.html
https://stackoverflow.com/questions/14492284/center-a-map-in-d3-given-a-geojson-object
Conversion from Shapefile to TopoJSON format results in incorrect map when rendered using d3js
http://milkator.wordpress.com/2013/02/25/making-a-map-of-germany-with-topojson/
and Mike's original example
http://bost.ocks.org/mike/map/
Maybe it is because i am not using EPSG:4326, but this example says that you do not have to use a projection
http://bl.ocks.org/mbostock/5557726
The data can be found here - https://www.dropbox.com/s/dfxxi4xwrsdfa8t/southeast.json
And this is my code
<script src="http://d3js.org/d3.v3.min.js"></script>
<script src="http://d3js.org/topojson.v1.min.js"></script>
<script>

var width = 500,
    height = 500;

//null projection as my ESPG is not supported
var path = d3.geo.path()
    .projection(null);

var svg = d3.select("body").append("svg")
    .attr("width", width)
    .attr("height", height);

svg.append("rect").attr('width', width).attr('height', height)
    .style('stroke', 'black').style('fill', 'none');

//d3.json("southeast.json", function(error, south) {
  //console.log(south);
//});

d3.json("southeast.json", function(error, uk) {
  svg.append("path")
      .datum(topojson.feature(uk, uk.objects.southeast))
      .attr("d", d3.geo.path().projection(null));
});

I first tested with sending the data to the console and that worked fine and told me that the object is "southeast".
If I try this I get an n is not defined error in Firebug.
I then tried
// d3.json("southeast.json", function(error, south) {
//   svg.append("path")
//       .datum(topojson.feature(south, south.objects.southeast))
//       .attr("d", d3.geo.path().projection(null))
//       .style("fill", "red")
//       .style("stroke-width", "1")
//       .style("stroke", "black");
// });

Which seems to work as there are no Firebug errors, but no data is visible
Then I tried
// d3.json("southeast.json", function(error, uk) {
//   svg.selectAll(".southeast")
//       .data(topojson.feature(uk, uk.objects.southeast).features)
//     .enter().append("path")
//       .attr("d", path)
// });

And the same thing happens
With no errors I am struggling to understand what I have done wrong
If I could just make this example work in EPSG:27700 with my topojson it would be amazing
http://bl.ocks.org/mbertrand/5218300
Thanks
EDIT
One of the other answers was to use a simply getJSON function and add the data as geojson in OpenLayers. I have just tested it using the same topojson:
$.getJSON('southeast.json', function (data) {
            fusiegem = topojson.feature(data, data.objects.southeast);
            var geojson_format = new OpenLayers.Format.GeoJSON();
            var vector_layer = new OpenLayers.Layer.Vector();
            osMap.addLayer(vector_layer);
            vector_layer.addFeatures(geojson_format.read(fusiegem));
        });

And that successfully loaded, so it seems there is nothing wrong with my topojson file, it is the D3 bit that does not seem to work correctly


texts: Anyone know if a package for processing SAR data in R exists? Such as doing calibration, terrain correction, classification on e.g. Sentinel 1 data?


texts: I have downloaded QGIS version Lisboa 1.8.0. I created a server-client architecture where one PC acts as server, in which QGIS is installed, while the others act as clients, which don't have a CPU. I am able to run MS Office, videos, etc. from the client PCs, but QGIS runs only on the server. What could be the possible problem and solution for it..? Is there any problem of user rights or codes..??


texts: I am relatively new to qgis. I have GPS coordinates of two types of sites in csv format, (water point site and sanitation site). I wanted to find out if it is possible to create a 100 meter radius on all the water point sites to find out how many of the sanitation sites are outside the 100 meters.


texts: Identifying Compression Ratio in a ECW raster file

texts: I'm trying to identify and verify the compression ratio of an ecw raster file and I cannot find the information under properties > source tab. I'm trying to using ArcMap 10.5 to do so.


texts: arcpy interate row, in a python script? SQL Select expression not working

texts:
I am creating a python script where I am trying to iterate over point feature class. I want to create seperate Feature classes from each individual day. So select daynumber and then export it to a new FC with an unique daynumber as name.

I was unable to find my answer online. So I decided to ask my question here.

ArcGIS 10.2

My script looks as follows:

test1 = "D:\\Mijn_Doc\\Test_Geodatabase.gd\\test1"
dayno = 2618
with test1 as rows:
    for row in rows:
        query = '"daynumber = '+ str(dayno) + '"'
        arcpy.Select_analysis(test1, 'D:\\Mijn_Doc\\Test_Geodatabase.gdb\\select'+str(dayno), query)
    dayno +=1

Gives the following error:

  arcgisscripting.ExecuteError: ERROR 000358: Invalid expression "daynumber = 2618" Failed to execute (Select)

The selection looks fine, If I build it in Model Builder and export it to a Python script the expression looks the same as it looks in this Error message.


texts: Using QGIS API and Python, to return latitude and longitude of point?

texts:
I have a point layer, that I'd like to return the Longitude, latitude using the QGIS API and Python.

How do I use QgsPoint to do that?


texts: Legend label failed to send to print service using ArcGIS javascript?

texts:
I am try to print my map with legend to PDF format; however, the legend label in the PDF is saying Override 1 instead of the label I specified in class breaks.

There are other maps I can print successfully with the correct labels used and the only difference I can find between them is that the ClassBreaksRenderer constructor is called differently. While the maps that can successfully print with the legend label is passing string to the attributeField parameter likevar renderer = new ClassBreaksRenderer(null, 'break_column_name');the one failed to print the legend label is passing a function to the parameter likevar renderer = new ClassBreaksRenderer(null, calculate_break);

I have searched the web and as someone suggested, looking at the JSON that sent to the print service, I found that that no information about the breaks (i.e. classBreakInfos) are sent to the service.

Does anyone know why it is the case and how can I fix it?

I am using the 3.10 API.


texts: Set of standard accessibility guidelines for online maps?

texts:
Recently I've been designing maps without navigation tools like the Zoom In/Out and Pan tools, or arrows to pan north, pan west, etc.

Instead, most modern mapping APIs support the direct use of the map itself to handle navigation - drag the map to pan, double-click to zoom in, scroll the mouse wheel, shift-drag to zoom, etc.

Is there a set of standard accessibility guidelines for online maps?

Are there valid reasons for needing a "click to pan" control, or any of the other old-school navigation controls? For example, are they used by screen readers or other accessibility aids?

I note that Google Maps has a Pan control, even though it's not strictly necessary.

Thanks.


texts: Adding hole to polygon using ArcGIS Online Web Mapping Application

texts:
Is it possible to add a hole / donut / doughnut to a polygon when using an ArcGIS Online Web Mapping Application?
I have found the Edit widget but whilst is has a Cut, Reshape and Union function, it does not appear to allow cutting a hole in a polygon.


texts:
I am trying to build a tool that will precisely (10digit, eg: AA1234512345) identify British National Grid (BNG) by clicking on a map. The ArcGis default tool can only display 6digit (eg. AA123123).
I have already written a code that calculates required BNG for point/line/polygon within a shapefile/feature class and recording it in an attribute table (I hope to cannibalise it).

I know that Arc 10.1 introduced something like onMouseDownMap [http://support.esri.com/es/knowledgebase/techarticles/detail/40730] (and Python Add-in Wizard) which seems to solve my problem but I am using Arc 10.0…

Is there anything in Arc 10.0 similar to onMouseDownMap? How to read coordinates (X, Y) from mouse click.

In brief my tool should work like:

Click on ‘identify’ button --> click on a map --> pop up window with BNG

My knowledge of python is rather limited… but I am fighting.


texts: How does QGIS handle opening files and how can I improve this?

texts:
We are using shapefiles with a large attached database to it (162 Mb for 133000 features). Opening this shapefile on a remote server takes about 5 minutes whereas DivaGIS takes 7 seconds for the same file.

My suspicion is that QGIS loads the entire shapefile to client and stores it in a temporary directory before rendering it whereas DivaGIS picks up the geometry, renders it and opens the .dbf file on the remote server.

My question is: "How does QGIS handle files and how can I improve this and / or can it be improved?"

(I have once read an article "QGIS on steroids" but I have not seen any improvement since QGIS 1.6)

Working with Postgis improves this a little bit but it doesn't come anywhere near the speed of DivaGIS.


texts: Setting up travel direction when creating network dataset?

texts:
I want to conduct service area analysis with a roadway network to create travel time contour map. The roadway netwotk is a linear shapefile which stores travel direction (N/S/W/E) and travel time for each roadway segments. So for a two-way road, there will be two lines with opposite travel direction and different travel time. So different lines would be used in calculating travel time from and towards the facility.

However, I don't know how to setup travel direction when creating the network dataset by using the existing travel direction.

Maybe I can set the travel direction of all lines as one-way and let GIS to determine the travel direction? However, the digitized direction of the lines are not the same with the travel direction.


texts: MakeRasterLayer failing with error codes 000622 and 000628?

texts:
I am importing a raster file and then saving the layer file. This script has worked fine previously, but after changing some input file paths, it has stopped working. Arc finds the inputs, but then has an error code at the make raster layer point related to the input parameters: error codes 000622 and 000628. After making the raster layer I save the layer file, and do other steps, but I have only included the script to the point where it fails.

    # Import arcpy module
    import arcpy

    # Import system modules
    import sys, string, os, arcgisscripting

    # Create the Geoprocessor object
    gp = arcgisscripting.create()
    gp.Overwriteoutput = 1

    # Set Workspace
    arcpy.env.workspace = r"S:\Map_Creation"

    # Set overwrite option
    arcpy.env.overwriteOutput = True

    # Local variables:
    InputRaster = "C:\\Asc_Conversion\\rasters\\r001d1b1_dr0"
    ExtentFile = "S:\\01.RIM\\Final_Extent_Outputs\\r001d1b1_ext_r0.shp"
    RasLayerTemp = "Depth"

    # Process: Make Raster Layer
    arcpy.MakeRasterLayer_management(InputRaster, RasLayerTemp, "", ExtentFile, "")


texts: How to create a point along a line given distance

texts:
I have been given a non-spatial database which contains the names of lines which I have a shapefile for. The database also has names of points and the distance of each point from the origin along the route of these lines.
How would I create a point shapefile which locates each of the points?

Of course I would like to use QGIS, or GRASS and as a last resort I can get access to ArcGIS. I expect I will have to carry out a step to specify which end of the line is the origin.


texts: How to get WKT from SHAPE column in the ArcSDE geodatabase?

texts:
I have an ArcSDE geodatabase running on top of SQL Server 2008. A table for a feature layer was added from ArcCatalog. The SHAPE column is a type of int which was automatically converted from ArcCatalog.

My question is:

In SQL Server Management Studio, what is the best SQL syntax to output WKT for each of the records in the table?

p.s. [SHAPE].STAsText() doesn't work in

SELECT [OBJECTID] ,[SHAPE], [SHAPE].STAsText() FROM [sde].[sdl_admin].[POINT_LOCATION]

Cannot find either column "SHAPE" or the user-defined function or aggregate
SHAPE.STAsText", or the name is ambiguous.


texts: -torque-resolution equivalent for static CartoDB maps

texts: I'm trying to create a Torque map with tooltips. I couldn't see a way to create tooltips on Torque layers, so as a workaround I've created a separate static layer than can be toggled. Toggling this layer sets the SQL to query the data for a point in time corresponding to the current Torque time. The Torque layer does some clustering/aggregation based on the -torque-resolution CartoCSS parameter. The problem I am now encountering is that the layers looks slightly different because the static layer doesn't do this clustering. Is there any way to replicate the behaviour of -torque-resolution CartoCSS on the static layer with SQL queries or otherwise? I actually want this clustering, so setting this to -torque-resolution: 1 won't help me in this case.


texts: Slow Labeling in ArcMap 10 (ArcView License)

texts:
We recently updated most of our users to ArcGIS10. One of my users is having an issue with incredibly slow labeling in ArcMap. The shapefile geometry displays right away, but the labels take minutes to show up at all. The shapefile is of all the parcels in a county and it seems like maybe ArcMap is labeling all the parcels, not just the ones in the current extent. It works just fine in ArcMap with an ArcInfo license. I don't have access to an ArcView 10 license to troubleshoot myself.

Are there any significant changes between ArcMap 9 and 10 for ArcView licenses that might cause this and is there's any way to fix it?


texts: Loading GeoJSON by AJAX method in OpenLayers, centering the map automatically

texts:
I'm trying to center the map on the geometries that are on the map. To encompass all but until now I have not succeeded in how ...

It centers further than geometries :(

  if (ol.proj.proj4 && ol.proj.proj4.register) { ol.proj.proj4.register(proj4); }

var features;

$.ajax({
    type: "GET",
    url: "get_geojson.php",
    dataType: 'json',
    async: false,
    success: function (response) {
       features = new ol.format.GeoJSON().readFeatures(response);
    }
});

  features.forEach(function(feature){
    feature.setId(undefined);
    feature.getGeometry().transform('EPSG:31700', 'EPSG:3857');
  });

  var map = new ol.Map({
    layers: [
        new ol.layer.Tile({
            source: new ol.source.OSM(),
        }),
        new ol.layer.Vector({
            source: new ol.source.Vector({
                features: features
            })
        })
    ],
    target: 'map',
    controls: ol.control.defaults({
            attributionOptions:
                ({
                    collapsible: false
                })
        }).extend([
            new ol.control.OverviewMap(),
            new ol.control.FullScreen(),
            new ol.control.ScaleLine()

        ]),
        interactions: [
            new ol.interaction.DragRotateAndZoom(),
            new ol.interaction.MouseWheelZoom(),
            new ol.interaction.DragPan(),
            new ol.interaction.DoubleClickZoom()
        ],
    view: new ol.View({
        center: ol.proj.transform([698282.236, 377375.758], 'EPSG:31700', 'EPSG:3857'),
        zoom: 14,
    })
  });


texts: Are there any shapefile-based GPS data recorders?

texts:
I am looking for a data recorder that can be programmed to send a serial data request to a field computer and then accept and record an incoming data string and tie it to GPS location and then save the data as a shapefile at the end of the logging period.  Any recommendations or places I can start looking?

I tried to keep my original post relitivley general to reduce confusion. The field computer I refered to in my original post is a spray control console that controls spray application in agricultural fields.  it has a speed sensor, flow meter, control valve, and gps receiver attached to it.  its purpose is to maintain a given application rate regardless of changes in speed.  The console has in the past had an onboard CFcard reader/data recorder that initiated a shapefile on power up of the console and closed the file on power down of the console. on the next power up a new shapefile. The CF card readers are no longer being produced so the consoles are now being produced with out them, thus creating the need to come up with an alternative.  The console has a serial comunications connector on it that will output a string of data.  the current onboard recorder did this at programamble intervals, internally and also collected the location information from the connected gps receiver and ultimately created a shapefile consisting of a SHP, DBF, SHX files.  As this is what we are accustomed to working with, and my customers are fairly tech handicapped, it seemed easiest to keep things as similar as possible.  Once the files have been produced they open the files with a GIS software package to view the application data, usually in a color coded map format, but also open the DBF file up directly to view the database lines directly.
     So what ideally I need is a recorder with 2 serial connectors.  one to receive the console data string and one to accept the GPS string,  on completion of the job the recorder will need to save the collected data to a Shapefile.
     If there are simple way to collect the data from the unit and location stamp it in a way that can be viewed by a shapefile viewer that would be fine to.
Thanks


texts:
I use Proj4j for a WebGIS and I have to create a new ESPG.
I tried on QGIS to convert coordinates from a local system (in Greece) to a world coordinates system.

First, I have to do a 2D transformation: 2 translations (Tx ans Ty) and a rotation (Rz) of about 0.58931 radians around the "z axis" (in the plane of the local system) between the local coordinates and the Greek coordinates (EGSA87).

So I think that the rotation matrix is: [cos(Rz),-sin(Rz);sin(Rz),cos(Rz)] and without the cos and sin I have [0.8313242,-0.5557878;0.5557878,0.8313242].

Then, I can reuse the parameters of the EPSG:2100 definition:

+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=500000 +y_0=0 +ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs

I think I can insert the translations (Tx and Ty) in the x_0 and y_0, but I don’t know how to insert the rotation matrix. How can I insert a transformation with a matrix in a new EPSG?

I found on the internet that the +xform= command seems good, but I don’t know how to use it. I tried:

+proj=tmerc +lat_0=0 +lon_0=24 +k=0.9996 +x_0=-869388.894 +y_0=-4260031.458 +xform=[0.8313242,-0.5557878,0.5557878,0.8313242] +ellps=GRS80 +towgs84=-199.87,74.79,246.62,0,0,0,0 +units=m +no_defs

But apparently it isn’t a good solution. So, does anyone know how I can use this command in my new EPSG or where I can find the proj4 documentation?


texts: Converting Esri Leaflet featureCollection to JSON?

texts:
I've followed the Esri Leaflet sample code found here:

http://esri.github.io/esri-leaflet/api-reference/tasks/query.html

//Querying features near a latlng
var latlng = L.latLng(45.51, -122.70);

var query = L.esri.Tasks.query({
            url:'http://services.arcgis.com/rOo16HdIMeOBI4Mb/arcgis/rest/services/stops/FeatureServer/0'
});

query.nearby(latlng, 500);

query.run(function(error, featureCollection, response){
console.log('Found ' + featureCollection.features.length + ' features');
});

So now that I have the features any Idea of how to export them out as JSON or GeoJSON? I need them as input for another app.  So far the only thing I found you could do is count them.


texts: Meaning of empty quotes in select by attribute query?

texts:
Can someone please help explain if the empty '' means something?

I am wondering if it might be inserted by typo or if you use these empty quotes to locate empty fields within the attributes.


texts:
I'm converting PostGIS data to PDF via GDAL 2.0.2 with below.

ogr2ogr -f "PDF" mypdf.pdf PG:"host=domain.com user=me password=mypassword dbname=mydb"  \
 -sql "SELECT geom,  state_name,  sub_region, state_abbr FROM  states"

It outputs OK, but I can't seem to find the method to get the output PDF to display the state_name column in the tree and not just feature1, feature2, etc.

Working now:

ogr2ogr -dsco OGR_DISPLAY_LAYER_NAMES="States" \
        -dsco OGR_DISPLAY_FIELD=state_name \
        -f "PDF" \
        mypdf.pdf \
        PG:"host=domain.com user=me password=mypassword dbname=mydb" \
        -sql "SELECT geom, state_name, sub_region, state_abbr FROM states"


texts: Compute geometric mean of a field (statistical computation)

texts:
I am trying to use the geometric mean instead of the arithmetic mean to aggregate a field in ArcGIS-online. Arcade does not offer a function for the geometric mean (like Average($feature["Result"] for the arithmetic mean), so I try to write it myself:
var product = 1;
var counter = 0;
for(var k in $feature["Result"]) {
    Console('looping!'); // not printing in the Messages tab
    product *= $feature["Result"][k];
    counter += 1;
}
var power = 1 / counter;
Console(counter); // printing in the Messages tab
// the counter is never updated, meaning the loop does not run...
// resulting in a division by 0.
var geomean = Pow(product, power);
return geomean;

But it does not work. I realize that the fundamental problem here is: How to loop through the value of one field? I do not want to loop through several field names, I would like to loop through all the values of this "Result" field. Is it possible? There is a fundamental understanding of how this works that I cannot grasp...
Thanks for your help!


texts: How can I display this polygon in Google Maps, properly

texts:
NOTE: this question is sorta like a follow-on question to this original GIS Question I asked. In that question, the main suggestion was that I was trying to show two polys. Now, I'm asking about a SINGLE poly.

I have the following GEOGRPAHY of Los Angeles in my Sql Server 2008 :-

Sql script to repo it (or in RAW sql format)

It's bascially ..

-- Shape of Los Angeles.

DECLARE @WKT VARCHAR(MAX) = 'POLYGON ((-118.66815700000011 34.181235000000029, ... ))
SELECT GEOGRAPHY::STGeomFromText(@WKT, 4326)

Now when I just try and render this out, It now looks like this (NOTE: this is different, but still not correct .. compared to my previous GIS Question).

Now, in the previous question people were suggesting that the previous data was a MULTIPOLYGON so there was some bad stuff I must have been doing, trying to render out holes in polys.

This time, this is a single polygon. So i'm assuming the hole theory still doesn't work. Now, I understand that by looking at this image, there's a number of holes ... but the Well Known Text data which Sql Server gives me is a POLYGON which to me means .. one single shape. So i'm not sure how I can fix this :(

This is why i'm asking the a similar question to my previous one.

Can anyone help? Can they see the data I've posted and see what happens in other systems? (I don't have any other software besides Sql Server 2008).

Maybe I need to extract all my data from Sql Server 2008 .. and clean it up (whatever that means) through another system. For example. If Sql Server is saying that the shape of Los Angeles is a single Polygon .. but when cleaned up .. it's actually 7 Polygons .. and then those shapes then get rendered correctly .. then we have a solution!

Any suggestions?

Can -anyone- render this shape onto a Google map?


texts:
I am a complete novice.

I have the following .prj file:

PROJCS["EO_Lambert_Conformal_Conic",GEOGCS["GCS_North_American_1983",DATUM["D_North_American_1983",SPHEROID["GRS_1980",6378137.0,298.257222101]],PRIMEM["Greenwich",0.0],UNIT["Degree",0.0174532925199433]],PROJECTION["Lambert_Conformal_Conic"],PARAMETER["False_Easting",1000000.0],PARAMETER["False_Northing",0.0],PARAMETER["Central_Meridian",-84.0],PARAMETER["Standard_Parallel_1",44.5],PARAMETER["Standard_Parallel_2",54.5],PARAMETER["Latitude_Of_Origin",0.0],UNIT["Meter",1.0]]

I need to get the EPSG code for this, but when I try http://prj2epsg.org/search I keep getting the following error

  Recursive call while creating a 'HTTP_AuthorityFactory' object.

I have searched all over, but cannot figure this out.  I have tried http://spatialreference.org/ and found http://spatialreference.org/ref/sr-org/7149/ which I think is the same .prj file, but 7149 is not a valid SRID according to PostGIS.


texts: Converting GeoJSON to shapefile in Python with Fiona: problem with Geometry Collection

texts:
I have a function, using the Fiona package, to convert a GeoJSON to their respective shapefiles (going to figure out later to get merge them back together later on) However the problem is the GeoJSON, multiple geometries. And this causes the following error. So according to this I need to create a layer for each:
Geometry type of 'Geometry Collection' not supported in shapefiles.  Type can be overridden with a layer creation option of SHPT=POINT/ARC/POLYGON/MULTIPOINT/POINTZ/ARCZ/POLYGONZ/MULTIPOINTZ/MULTIPATCH.

What should I do here?
Also of note this GeoJSON can only have two types (technically three with geometry collection) Polygon & MultiPolygon. So everything should turn up in one shapefile later on.
import fiona
from fiona.crs import from_epsg

 def convert_json_to_shapefile(fc_file):
       total_start_time = start_time_measure()
       start_time = start_time_measure(
       'opening geojson file')
       export_folder = os.path.join("exported_data","geographic_information","details","shapefile")
        # First, we load the JSON file, and create a var to hold the data
       with open(fc_file) as f:
            fc_data = json.load(f)

        # Create export folder if it doesn't exist
       if not os.path.exists(export_folder):
            os.makedirs(export_folder)

        # Iterate over each feature
       total_num_features = len(fc_data['features']) - 1
       for num, feature in enumerate(fc_data['features']):
            print(f"Parsing feature {num} of {total_num_features}")
            # Let's get some data from each feature
            gipod_id = feature['properties']['gipodId']
            reference = feature['properties']['reference']  # Replace slashes with dashes for path
            geo_type = feature['geometry']['type']
            geometry = feature['geometry']

            # Split the EPSG string, get the last value, convert to int and generate a Fiona EPSG
            split_epsg = int(feature['geometry']['crs']['properties']['name'].split('::')[-1])
            feat_epsg = from_epsg(split_epsg)

            # Print out the data
            print(f"GIPOD ID: {gipod_id} - Ref: {reference}")
            print(f"Geometry type: {geo_type}")
            print(f"Spatial reference: {feat_epsg}")
            print(feature['properties'])
            # Create a geometry schema
            geo_schema = {
                'geometry': geo_type,
            }
            # Export filename. For now we'll just use the gipod ID and ref
            print(str(reference))
            export_file = os.path.join(export_folder,
                                       f"Shapefile_ID{str(gipod_id)}.shp")
            print(f'Writing {export_file}...')
            with fiona.open(export_file, mode='w', driver='ESRI Shapefile', schema=geo_schema, crs=feat_epsg) as output:
                output.write(geometry)
            print('----------')

   print(f"Done! You can find your shapefiles in the {os.path.abspath(export_folder)} folder. :-) ")


texts:
I am trying to import the pymodis library within python and I keep getting error that

"GDAL installation has no support for HDF4, please update GDAL"

I am running the latest version of GDAL [Version: 1.11.2+dfsg-1~exp2] and tried runnign the update command too. But still the same error persists.

I am running Ubuntu 15.04 and using the terminal to run  python.


texts:
I am trying to delete a domain I am using from ArcCatalog -> Feature Class Properties-> Subtypes-> Domains. But when I try to do that I get this error message:
"Failed to delete a domain from the database.
The domain is used by an attribute rule."

I've checked if there are any locks and there aren't.

Does anyone know  why I can't delete it?

Thank you

I use ArcGIS 10.1. Standard on an ArcSDE Postgres Geodatabase.

Update: the problem degenerated a lot :) I was able to delete the field corresponding to the domain, but I still cannot delete the domain. Same error. Plus, when I want to disable the existing subtype I get this message  "Failed to change the subtyoe field. LInestring or poly boundary is self-intersecting". I assume it's just one single problem that messed up the subtype and the domains and I get different error messages.


texts:
I created an ArcSDE geodatabase via create enterprise database tool. I used MSSQL Server. But normally there has to be an icon which is exactly the same one left side of  "SLK sql to sa.sde" text under the layers. You know, underlined of database icon.  Because in the past, when I dropped the feature class, there was that icon which can be selected version user type ( Default, child etc. ) but know there is only a layer and just like ordinary feature class. I want to select version ( child /parent ) but there is no selection.

What is my mistake?


texts: Subsetting one raster based on attributes from another

texts:
I am trying to expand existing 'tidal complexes' to include the adjacent 100-year floodplain areas.  Both are rasters, and so far I've used Euclidean Allocation to identify the tidal complex closest to each pixel in the study area and extracted the output using the floodplain as a mask so that I have a raster of floodplain area labeled by the closest tidal complex.

However, I only want to include floodplain areas contiguous to the tidal complexes - some of them are farther inland and not contiguous.

I region grouped the floodplain area labeled by tidal complex output so that each contiguous area has a unique value, and then extracted that output using the original tidal complexes as a mask.  The resulting raster has the unique identifiers for the contiguous floodplain areas underlying each tidal complex.

I'm stuck at this point; I need to somehow extract the floodplain areas contiguous to the tidal complexes (my current raster's extent is just the tidal complexes themselves).  I originally thought that an inner join between the region-grouped floodplains and the tidal complexes identified by region-grouped floodplain ID would do the trick, but it seems to alter just the raster attribute table, not the raster itself - the attribute table only includes the region-grouped floodplains contiguous to the tidal complexes, but the raster has all of the region-grouped floodplains.

My other idea is to select the region-grouped floodplains by attribute, based on the values in the tidal complexes by region-grouped floodplain ID, but there are >200 values, and I don't know of an efficient way to extract them to use in a query.

I am working in ArcGIS Pro ModelBuilder so far, but could switch to python if necessary.


texts: Download Copernicus Marine variables with Package RCMEMS - Error code 9

texts:
I'm using the source code RCMEMS (https://github.com/markpayneatwork/RCMEMS) to download several marines variables and I have a error despite having respect syntax package.
To biginning, according to my python terminal I use Python 3.6.0 and motuclient-python v1.8.6:
F:\PNB_ENTROPIE>python --version

Python 3.6.0

F:\PNB_ENTROPIE>python -m motuclient --version

motuclient-python v1.8.6

For example, I tried to download SST variable (https://resources.marine.copernicus.eu/?option=com_csw&view=details&product_id=SST_GLO_SST_L3S_NRT_OBSERVATIONS_010_010) with VScode and R language :
devtools::install_github("markpayneatwork/RCMEMS", force = TRUE)

cfg_manuel <- RCMEMS::CMEMS.config(
                    motu = "http://nrt.cmems-du.eu/motu-web/Motu",
                    python = "C:/Users/romai/AppData/Local/Programs/Python/Python36/python.exe",
                    script = "C:/Users/romai/AppData/Local/Programs/Python/Python36/Scripts/motuclient-script.py",
                    user = "rfernandez1",
                    auth.mode = "basic",
                    pwd = "XXXXXXX",
                    service.id = "SST_GLO_SST_L3S_NRT_OBSERVATIONS_010_010-TDS",
                    product.id = "IFREMER-GLOB-SST-L3-NRT-OBS_FULL_TIME_SERIE",
                    longitude.min = "10",
                    longitude.max = "155",
                    latitude.min = "-55 ",
                    latitude.max = "31",
                    date.min = "2020-12-16",
                    date.max = "2020-12-31",
                    variable = "sea_surface_temperature",
                    out.dir = "F:/PNB_ENTROPIE/Variables",
                    out.name = "Test_RCMEMS"
)

When I want download my variable stocked in cfg_manuel  object :
RCMEMS::CMEMS.download(cfg_manuel)

I have this output:

2021-01-19 15:00:32.157 [ INFO] Asynchronous mode set
2021-01-19 15:00:32.160 [ INFO] Requesting file to download (this can take a while)...
2021-01-19 15:00:34.557 [ERROR] Execution failed: junk after document element: line 30, column 1
2021-01-19 15:00:34.558 [ INFO]  . code  9:
Error in (function (x, out.dir = NULL, out.name = NULL, date.min = NULL,  :
  Error in running CMEMS download command.

What is this error and where does it come from?


texts: Enabling CORS in GeoServer 2.8.2 or how to use JSONP in GeoServer?

texts:
I want to identify number of layers in to workspace in GeoServer for that  I am using below REST API:

http://admin:geoserver@localhost:8080/geoserver/rest/workspaces/gisc/featuretypes.json?callback=?

My application is running on IIS so I am using JSONP. I have enabled JSONP in GeoServer but I can not get output, always getting error "**Uncaught SyntaxError: Unexpected token :"

I tried:

var url  = "http://admin:geoserver@localhost:8080/geoserver/rest/workspaces/gisc/featuretypes.json?callback=?";

$.ajax({
 type: 'GET',
 url: url,
 async: false,
    jsonpCallback: 'parseResponse',
    contentType: "application/json",
    dataType: 'jsonp',
    success: function(data)
    {
        console.log(data);
    },
    error: function (e)
    {
        console.log(e);
    }
 });

I have also tried another combination but it won't work I am always getting same error.

I have also tried to enable CORS in Jetty using "org\mortbay\servlets\CrossOriginFilter.class" it works with WMS and WFS but won't work with featuretypes.json.


texts: JOSM data to be exported and used in GIS

texts:
I need to use open street map data, specifically the turning lanes's data. I can very easily import it in JOSM but it would be useful for me if I had it o GIS format.

I am not very familiar with JOSM but it seems that I can't find a 'to export' turning lanes data, or another data from JOSM.

Is there any way to do so?


texts:
when I use ArcGIS I normally use CRS: Cylindrical Equal Area (World) WKID: 54034

In QGIS I cannot find the same projection. When I re-project my file in ArcGIS and open it in QGIS it says I am using a custom CRS:

Generated CRS (+proj=cea +lon_0=0 +lat_ts=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs)

Is it true that QGIS doesn't have this CRS? I will use my custom CRS for now but for future projects my main question is: What is a suitable equal area projection with global coverage in QGIS? I was thinking about Lambert but can't find the right flavor.


texts: Option for 360 degree movie embedded in Open Source GIS software?

texts:
Is there any Open Source GIS software available to embed 360 degree movies with the map?

For a COTS example, there is the Active Link Vision (ALV) extension for ArcGIS from IWANE.


texts: Sharing information about apps built with ArcGIS API for JavaScript?

texts:
I built an app with the ArcGIS Javascript API and made it available on GitHub.

Where would be an appropriate place to tell GIS users/developers about it?


texts:
Has anyone ever used the two of them together, and if so, any special caveat to be aware of ?

Furthermore, does it even make sense to use PostGIS / do TerramapServers have their own storage format instead ?


texts: find nearest point to a point in another table

texts:
This similar, but not quite the same as this question.

I need to pick, regardless of distance, the closest 'n' points to a given point in another table, using postGIS.

Imagine a table with city names, state names, and a centroid for the city itself, and another table with point information on various weather radar sites. Both tables re postGIS enabled.

PostGIS-2 has me confused; it's been awhile since I had to do something line this.


texts:
I am trying to build GeoDataFrame row by row, from an empty one. The equivalent using only pandas would be something like this :

df = pandas.DataFrame(columns=['a','b','c','d'])
df.loc['y'] = pandas.Series({'a':1, 'b':5, 'c':2, 'd':3})

(from [this answer])

My current method

So far, I build a Python list of dict with a specific structure, and then use it to create a GeoFataFrame, here is a complete example :

import geopandas as gpd
from shapely.geometry import  Point

my_dict = {
  '007': {
    'name': 'A',
    'lat': 48.843664,
    'lon': 2.302672,
    'type': 'small'
  },
  '008': {
    'name': 'B',
    'lat': 50.575813,
    'lon': 7.258148,
    'type': 'medium'
  },
  '010': {
    'name': 'C',
    'lat': 47.058420,
    'lon': 15.437464,
    'type': 'big'
  }
}

tmp_list = []
for item_key, item_value in my_dict.items() :
  tmp_list.append({
    'geometry' : Point(item_value['lon'], item_value['lat']),
    'id': item_key,
    'name': item_value ['name'],
    'type': item_value ['type']
   })
my_gdf = gpd.GeoDataFrame(tmp_list)
print(my_gdf.head())

Here is the result :

                    geometry   id name    type
0   POINT (2.30267 48.84366)  007    A   small
1   POINT (7.25815 50.57581)  008    B  medium
2  POINT (15.43746 47.05842)  010    C     big

What I am looking for

I would like to create an empty GeoDataFrame (my_gdf = gpd.GeoDataFrame()), and then fill it directly in the for loop, without using the temporary list after the loop (my_gdf = gpd.GeoDataFrame(tmp_list))

I think that a row by row building would, in this case, have better performance. It would also allow me to use the id key from my_dict as the GeoDataFrame index, so that the result would be :

                     geometry  name    type
007   POINT (2.30267 48.84366)    A   small
008   POINT (7.25815 50.57581)    B  medium
010  POINT (15.43746 47.05842)    C     big


texts: Investigating Avenza GeoTIFF Compression Algorithms

texts:
I was testing out different compression algorithms for import into Avenza and I found out that for the deflate algorithm, making the bit depth 24 decreased the file size and decreased the processing time.
However the confusing part was how to make the bit depth 24. Setting NBITS (in QGIS) to 24 made a 32 bit image (as per windows explorer right-click>properties). Yet setting NBITS to 25 or 17 made a 24 bit image.
The map was a general topographic map (without shading etc.)
My results were as follows

Packbits: processing = 34 seconds, size = 875MB
LZW:      processing = 41 seconds, size = 83MB
Deflate (32bit):  processing = 44 seconds, size = 64MB
Deflate 24bit: processing = 38 seconds, size = 60MB

The reason that setting NBITS to 17 created a 24 bit image is that: "Option name="NBITS" type="int" description="BITS for sub-byte files (1-7), sub-uint16 (9-15), sub-uint32 (17-31), or float32 (16)" /"
All the above was done in QGIS, meanwhile ArcGIS Desktop was able to make a 60MB LZW file that processed in 33 seconds so I guess QGIS's renderers aren't as optimised.
How do I reduce the file size as much as possible using QGIS for the importing of GeoTIFFs into Avenza?
EDIT: My eventual solution to this, as prompted by the ESRI support lady that called me up, was to jump onto ArcGIS Pro. From here I figured out that the solution is to:
Right click and zoom to your raster (this puts your view at 110%), then multiply the map scale by 0.909090909 (this brings your view in to 100%), then export the map as a deflate TIFF on max with tags with the resolution of the original raster (ArcGIS Pro seems to fix the weird pixel shifting that ArcMap had).
For QGIS it seems the best settings are Deflate at 24 bit (NBITS = 17) with some least significant bits removed (DISCARD_LSB = 4). Also always use NUM_THREADS = 8.
The result was 55mb in QGIS and 40mb in Arc


texts: Can't export GeoPDF on Mac from QGIS 3.12 (error "needs GDAL 3"), but GDAL 3 is installed already

texts: I am completely new to QGIS, and I am trying to make a custom Avenza map using it on Mac. When I follow the instructions (create print layout, export as PDF), I get all the way to exporting a PDF. The problem is that the Export GeoPDF checkbox is grayed out, and cannot be checked. There is a little note under the checkbox saying that I need GDAL 3 or greater. If I continue exporting, the PDF comes out without the GPS info in it. In the plugins section it shows that I have GDAL 3, but QGIS is obviously not using it. How do I fix this?


texts:
I am working with Linear Reference Systems and I know is very easy to get M coordinate (M-value) from a shapefile in ArcGIS. But, how to get the M-value in QGis?

If you have doubts about my issue, please look how ArcGIS works with M-value: http://resources.arcgis.com/en/help/main/10.1/index.html#//003900000006000000


texts: I have to use a qgis method (python, plugin, scripts) only to automate a shortest path calculation for multiple start and end points. How might i automate the shortest path road graph tool?


texts:
I'm trying to download QGIS 3.4. for mac. However web site says version 3.4., link name and downloaded version is still 3.2.x.

Anybody has to some issue? Do I anything wrong?


texts: Call a function multiple times with different argument for date

texts:
I need to call the same function 10 times, each time with a different year for the end date argument. For example, the end date needs to start as "2011-01-01", then "2012-01-01", then "2013-01-01", etc. How do I loop this function but change the one argument until the last end date is "2020-01-01"?
VI <- mt_subset(product = "MOD13Q1",
                band = "250m_16_days_NDVI",
                lon = 28.5826,
                lat = 1.4368,
                start = "2010-01-01",
                end = "2011-01-01",
                km_lr = 20,
                km_ab = 20,
                site_name = "Ituri test",
                internal = TRUE,
                progress = TRUE)


texts: Attribute values as input parameter in QGIS geopreocessing scripts

texts:
In QGIS geoprocessing scripts, I know how to connect an input parameter to a layer's attribute list. It would be like this:

##Input_layer=vector
##Input_field_name=field Input_layer

But, is there any way to connect an input parameter to an attribute's unique value list? Maybe something like this:

##Input_field_value=value Input_field_name


texts: ArcGIS vs QGIS vs PostGIS projection issue

texts:
Windows 10, ArcGIS Pro, QGIS 2.18, pgAdmin 4-postgreSQL

I have two (sample area) PostGIS tables in a PostgreSQL database:

building_polys (SRID: 27700) -OS MastermapTopo

address_points (SRID: 4258) -OS AddressBase

I have a project set up in British National Grid projection (27700) in ArcGIS and QGIS too. When my postgis layers are imported into ArcGIS I get the expected result where the point lies within the building polygon:

However when imported into QGIS, the result is wrong:

I believe this is due to the way QGIS re-projects layers from 4258 to 27700...

This in itself is not an issue, HOWEVER my problem is that using ST_WITHIN statements in PostgreSQL (using pgAdmin) seem to use the QGIS approach...

If I run:

SELECT t.*, s.uprn FROM public.building_polys t, public.address_points s WHERE ST_WITHIN(s.wkb_geometry,t.wkb_geometry)

It returns nothing. For comparison, if I add a 5 metre search distance on that (using ST_DWITHIN) I do get my expected results:

SELECT t.*, s.uprn FROM public.building_polys t, public.address_points s WHERE ST_DWITHIN(s.wkb_geometry,t.wkb_geometry,5)

This confirms that postgis recognises an offset just like QGIS.

I decided to apply a transformation of the point data, from 4258 to 27700, to see if that would shift the points into their correct position:

ALTER TABLE public.address_points
 ALTER COLUMN wkb_geometry TYPE geometry(Point, 27700)
   USING ST_Transform(ST_SetSRID(wkb_geometry, 4258), 27700);

This doesn't help though. It does transform the points to 27700 but they remain in the same position and more importantly, SQL statements still recognise the point as lying outside of the polygon:

How can I get my PostGIS PostgreSQL statements to recognise the spatial position of the data in the same way that ArcGIS does?


texts: How to determine which points fall within a certain value in a polygon in QGIS?

texts:
I have a shapefile that contains different levels for food security. On top of that, I have mapped points for certain office locations. This is shown below.

The country shapefile is extremely simplistic--the only information it contains are values for the food security index (values 1-5). Because of this, I am unable to perform any kind of spatial join.

I want to single out WHICH points exist in each different value, so I don't think point-in-polygon analysis is sufficient (as far as I can tell, this only gives me certain counts, not which points fall in each polygon).

What I've managed to do is select each food security value from the attribute table and then manually clip the point locations to each one. This is a bit clunky though, as I'm looking to do this for many countries. Is there a tool to achieve this more quickly?


texts: Accessing Spatial Analyst extension from stand alone script?

texts:
If I run this script:

import arcpy
arcpy.AddMessage(arcpy.CheckExtension("Spatial"))
from arcpy.sa import Con
from arcpy import env
fDir=r'd:\scratch\fdir'
outFolder=r'd:\aerials\images'
env.workspace = outFolder
fDir=arcpy.Raster(fDir)
Gter=Con(fDir> 1,2)
Gter.save("TEST")

from ArcGIS, the output is:

When I run it as stand-alone script, it is very different story:

How do I access Spatial Analyst extension from a stand alone script?


texts: Importing .csv files made from R:write.csv - some values are

texts:
When i make a data frame in R, e.g. like this:

x <- data.frame(sum=c(0.00, 0.00, 1363.31, 27663.58, 27663.58), ID=1:5)

I can save it with write.csv

write.csv(x, "x.csv")

But when i open the x.csv file in arcmap it reads the three last numbers in x$sum as < Null >. Why is this, and is there a way to ensure that ArcMap reads numbers correctly?

A workaround is to open it in excel and save it as a worksheet, but there has to be an easier way.

I dont think its the decimal separator - when I open it in notebook it looks like this:

"","sum","ID"
"1",0,1
"2",0,2
"3",1363.31,3
"4",27663.58,4
"5",27663.58,5


texts: arcpy.Intersect_analysis using Geometry object as input

texts:
Is there a way to run arcpy.Intersect_analysis using an arcpy.Geometry object as one of the input datasets?

I've tried different variations of the following, but I can't seem to get it to work.

with arcpy.da.SearchCursor(input_feature_class, ["SHAPE@"], "ObjectID = 100") as cursor:
    for row in cursor:
        geometry = row[0]

intersection = arcpy.Intersect_analysis([geometry,input_feature_class], arcpy.Geometry(), output_type="POINT")[0]

Results in:

 File "C:\Program Files (x86)\ArcGIS\Desktop10.3\ArcPy\arcpy\analysis.py", line 289, in Intersect
raise e
 ValueError: DescribeGeometry: Unexpected error

I'm already using a feature layer and SelectLayerByAttribute to achieve the same thing.  I'm curious as to why the above code doesn't work.


texts: ArcMap replica error number and description list

texts:
We have a replica between our sql server geoDB and another sql server geoDB. It is a one-way replica and our geoDB is the child. Recently an error occurs (error 000582). I searched for the cause of the error and got from ESRI forum which said:

  This error can occur if one or more of the feature classes and datasets are not versioned or don't have GlobalIDs."

Details check: https://geonet.esri.com/thread/52003.

  The feature classes and datasets are not versioned or don't have GlobalIDs,

May I confirm is the issue from the parent side?

Now the user want to see the official support for description of the error, but I failed to find it.

I want to check whether there is a source which list all the replica errors and the cause for them. I am using ArcSDE 10.1.


texts: Lining up three Layers that use same Coordinate Reference System (CRS) in QGIS?

texts: How do you make all borders of a layer equal? they should be equal but i don't know how, mCRS layer is all in philippine zone 5(EPSG:3125), with the same formula/equation, I'm new to QGIS and I'm using QGIS 2.4 Chugiak.


texts:
I want to work with a planar coordinate system, and as far as I know QGIS automatically puts me in a cartesian plan if I don't set up a CRS (Coordinate Reference System) for my project. That is exactly what I want.

The issue (may) begin when I import the raster image that I want to work with: it asks me to set a CRS for the image. The same happens when I create a vector layer for drawing features over my raster image: it asks me for a CRS for the vector too. The project remains without a CRS, but untill now I can not put any raster or vector without an CRS.

For tests purposes I set up EPSG:32722 - WGS 84 / UTM zone 22S as the CRS for both raster and vector. I used UTM because I know it is based on a cartesian plan (though I'm not sure it's the best to do).

My goal is to draw some polygons and calculate the area of that polygons. The scale is obtained from my raster (there's an object with known size in the image).

I can't measure correctly the size of my object without setting the CRS of the project to be the same of the raster and vector, so I set the CRS of the project from the layer.

It looks like the area is being correctly measured, but I can't know for sure if some transformation is going on.

That said, my questions are:

Is there a way to import raster or create vector layers without a CRS (am I missing something)?

Are my areas being calculated correctly this way?

Is alright to use UTM for this purpose or there is a better option (maybe define a custom CRS)?

(sorry for my bad English)

QGIS 2.2.0 64 bits on Win 7 64 bits


texts:
How to convert MapPoint to latitude/longitude in ArcGIS for Xamarine Android?
Below is my code.

MapPoint pointGeometry = new MapPoint(e.Location.X, Location.Y,SpatialReferences.Wgs84);


texts:
Is it possible to request the bounding box of a WMS from GeoServer?
In my case I have a WMS that I show with L.TileLayer.WMS in a Leaflet map. I like to center and zoom my map to the WMS.


texts: GML file converted to ESRI shapefile with ogr2ogr missing spatial reference information

texts:
I am trying to convert Sentinel-2 cloud masks from GML format to ESRI shapefiles. I have managed to use ogr2ogr in a python script as described in the answer to Using ogr2ogr to convert GML to shapefile in Python? to do this. I downloaded ogr2ogr from http://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/ogr2ogr.py and my Python call to that looks like this:

ogr2ogr.main(["", "-f", "ESRI Shapefile", 'output_shapefile.shp', 'cloud_mask.gml'])

However when I load the new shapefiles into ArcMap, it tells me that they are "missing spatial reference information". I've opened the GMLs in a text editor and I can see the EPSG code.

GDAL isn't creating a .prj file at all. If I try with '-lco SPATIAL_INDEX=YES' it creates a .qix file (which it doesn't otherwise), but the problem remains. The output of ogrinfo -al -so for the GML file is as follows:

  INFO: Open of 'cloud_mask.gml' using driver `GML' successful.
   Metadata:
    NAME=MSK_CLOUDS pixels mask from data-strip S2A_OPER_MSK_CLOUDS_SGS__20190403T
  205739_A019743_T13RGP_B00_MSIL1C

  Layer name: MaskFeature
   Geometry: Polygon
   Feature Count: 4
   Extent: (806760.000000, 3335340.000000) - (809760.000000, 3381900.000000)
   Layer SRS WKT: (unknown)
   gml_id: String (0.0) NOT NULL
   maskType: String (6.0)

Is there some way I can get ogr2ogr to project the shapefiles that it creates correctly, or another way I can project them afterwards using some extra Python code?